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PREFACE 


This report introduces aspects of the theory of continuous stochastic 
processes that seem to contribute to an understanding of turbulent disper- 
sion. It is expected that the reader has a knowledge of basic probability 
theory and some exposure to turbulence theory and the phenomena of turbulent 
transport. 

The report is addressed to researchers in the subject of turbulent 
transport, especially those with an interest in stochastic modeling. How- 
ever, the material presented herein deals with the theory and philosophy of 
modeling, rather than treating specific practical applications. The book by 
Csanady (ref. 1) is a good place to begin an exploration of applications, as 
well as to obtain a background to the contents of the present report. 
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CHAPTER I: OBJECTIVES AND INTRODUCTION 




The transport, or dispersion, of contaminant by turbulent convection is 
an archetypal continuous stochastic process. The theory of stochastic dif- 
ferential equations evolved, from early studies on Brownian motion, as a 
mathematical theory of continuous stochastic processes. It is natural to 
look to the theory of stochastic differential equations for insight into the 
phenomenon of turbulent transport. Indeed, some comprehension of the theory 
of continuous stochastic processes seems prerequisite to a basic understand- 
ing of turbulent dispersion. This report introduces aspects of stochastic 
differential equations that seem to bear on turbulent dispersion theory. 

The problem of turbulent dispersion can be posed as follows: 

Given a statistically prescribed field of turbulent velocity, describe 
statistically the evolution of a contaminant field from a known initial 
state. 

In this formulation of the problem the velocity field is given in whatever 
degree of minuteness is requested. From this information one wants to 
describe statistically the concentration field C(i<,t) of a contaminant, 
knowing it initially to be Cq(j(). This report is addressed entirely to the 

problem of determining either the mean value of this contaminant field 

"C(^,t) or the spatial moments of this mean contaminant field - the latter as 

functions of time alone. 

Two approaches can be taken to the analysis of turbulent dispersion, 
the Eulerian approach and the Lagrangian approach. In the Eulerian approach 
one formulates model conservation equations for fixed control volumes in the 
fluid; in the Lagrangian approach one formulates a random trajectory model 
of the motion of parcels of contaminated fluid. This latter approach will be 
followed herein. However, there is a close relation between the Eulerian and 
Lagrangian formulations: the Eulerian conservation equations determine the 

probability density function (pdf) of the random, Lagrangian, trajectories. 
For example, when the Eulerian model consists simply of an eddy diffusion 
equation (ref. 1, section 5.5), the Lagrangian model for which it determines 
the pdf represents trajectories by a continuous Markov process (chapter III). 
How could the evolution of mean concentration from an initial state be cal- 
culated from the statistics of random trajectories? 

Suppose one knows the initial distribution of contaminant to be Co(y) . 
Also, one knows an ensemble of random trajectories of fluid particles, each 
trajectory originating at some point yq . If it is assumed that fluid 
particles retain their initial concentration CQ(yo), then the mean concen- 
tration at the point (y,t) is just the probability that a parcel observed at 
position y at time t originated at some arbitrary point yo, times the 
concentration at that point Co(yo)» summed over all possible origins: 

■?^(y,t) = /* P(y,t; y^) CQ(yQ)dyQ 


(I-l) 



This equation uses the definition that the probability of a particle origi- 
nating at yo (actually, between yo - dyo/2 and yo + dyo/2) is 

P(y,t; yQ)dyQ (1-2) 

Note that in equation (I-l) the final position of the random trajectories is 
fixed at y, while the initial position yg is the random variable. (The 
random position yg would be calculated by following a trajectory from 
the known position y backwards in time to yg.) 

The derivation of equation (I-l) required that fluid particles conserve 
their initial contaminant concentration. To what extent is this valid? For 
any contaminant with finite molecular diffusivity k, fluid particles will 
not totally conserve their initial concentration. However, equation (I-l) 
describes only the average transport of contaminant; so this loss of contam- 
inant by molecular diffusion may not be important. Let us consider the rela- 
tive contributions of molecular and turbulent transport to mean dispersion. 

Molecular diffusion transports contaminant down concentration gradients. 
Hence, it acts primarily on small-scale features of the contaminant distri- 
bution, which have the largest gradients. These small-scale features are 
produced by an interaction between molecular diffusion and convection of 
contaminant by the small eddies of the turbulent velocity field. It follows 
that the length scale relevant to these small-scale features is the turbu- 
lence microscale (ref. 2, p. 67; the subscript < here indicating 
that the microscale is based on molecular diffusivity instead of on viscos- 
ity). Molecular diffusion, by smoothing these small-scale features, trans- 
ports contaminant a distance of order x< . Turbulent convection, on the 
other hand, transports contaminant a distance comparable to the turbulence 
integral scale (ref. 2, p. 20). The ratio of transport by molecular 
diffusion to that by turbulent dispersion is thus of order or of 

order (Pe)“l/2 (ref. 2, p. 68). The turbulent Peclet number is defined by 


Pe = u'a/< 


(1-3) 


where u' is the variance of the turbulent velocity. Hence, it can be con- 
cluded that equation (I-l) is formally valid when (Pe)“^'^ << 1. This 
report is concerned with high-Peclet-number and high-Reynolds-number turbu- 
lence, so equation (I-l) is justified. 

The problem of mean dispersion has been reduced by equation (I-l) to 
that of determining P(y,t;yg), or equivalently, to modeling the random tra- 
jectories for which P(y,t;yg) is the pdf. This can be done by representing 
the dynamics of fluid particles in turbulent flow by simple stochastic pro- 
cesses. The criterion such models must satisfy is that they reproduce cer- 
tain statistical features of the turbulent velocity field. As the problem 
of turbulent dispersion was posed above, the statistics of the velocity field 
are considered to be known in their entirety. However, the statistics most 
crucial to turbulent transport are the macroscopic length and time scales 


and the variance of the turbulent velocity and, in the case of nonhomogeneous 
turbulence, their spatial distribution. To make progress, models are based 
on only these restricted properties of the turbulence; for example, the eddy 
diffusion model requires knowledge only of the product u'x- of the length 
scale and velocity variance, which has the dimensions of diffusivity. The 
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present report is a discussion of concepts pertinent to stochastic modeling 
of turbulent dispersion. 

This report is not meant to be comprehensive, either in its treatment 
of stochastic differential equations or in its description of turbulent dis- 
persion. Stochastic differential equations are introduced in reference 3 
and treated more comprehensively in reference 4; for accounts of turbulent 
dispersion see references 1, 2, and 5. In chapter II a very nonrigorous 
review of the theory of stochastic differential equations is provided, touch- 
ing on those aspects of this theory that seem relevant to understanding and 
analyzing turbulent dispersion. Chapters III and IV deal specifically with 
models for the random trajectories of fluid particles through turbulent flow. 
If these trajectories are considered at times that are long compared with 
the macroscopic time scale Tl of the turbulence, they can be modeled by 
continuous Markov processes, or as they are alternatively called, by diffu- 
sion processes. Chapter III discusses eddy diffusion models; chapter IV 
discusses dispersion at times comparable to 1 \_ and uses the Langevin 
equation for a random trajectory model. 

To avoid later confusion, this introduction ends with a discussion of 
some notation to be used in this report. When a letter such as y is used 
as a dependent variable to denote a random trajectory, it will be written as 
a capital Y(t), with or without explicit indication of time dependence. 

When it is used as an independent variable, denoting a spatial coordinate, 
it will be written as a lower case y. The average of a random variable is 
signified by an overbar: Y(t) is the mean trajectory of Y(t). Of course, 

averaging can be carried out by integrating over a probability density if 
necessary: 


-f. 


GO 


y P(y,t)dy 


(1-4) 


This illustrates the use of upper- and lower-case letters. The function 
P(y,t) is the pdf of observing a particle at various points y at a certain 
t; y and t vary independently here. However, if at each time the m^an 
position of particles is computed, the result is the mean trajectory Y(t). 

In chapter II a time-dependent stochastic process, called the Wiener 
process, is introduced and denoted by W-^. By convention the time depen- 
dence here is indicated with a subscript: W-j- is synonymous with W(t). 

It might be added that if the randomness in Y is due to a dependence on 
the Wiener process, the average value of Y could formally be computed as 



Y(w^,t)P(w^,t)dw^ 


(1-5) 


In equation 
variable of 
to indicate 


(1-5) the lower-case w^ is used to indicate the independent 
integration; the subscript is for identification with W^, not 
time dependence. 
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CHAPTER II: REVIEW OF STOCHASTIC DIFFERENTIAL EQUATIONS 


The theory of stochastic differential equations was developed in the 
first half of this century by physicists and mathematicians. The physicists 
trace the origin of their research to Einstein's work on Brownian motion and 
include Langevin, Fokker, Planck, Smoluchowski , Ornstein, and Uhlenbeck 
(refs. 6 and 7). This line of research had as its impetus the use of con- 
tinuous stochastic processes to model random physical systems. The mathema- 
ticians' line of approach was to investigate formally the nature of continu- 
ous stochastic processes, including their relation to properties of partial 
differential equations. This body of work originated with Bachelier (who 
incidentally was investigating the fluctuation of stockmarket prices) and 
includes contributions by Wiener, Kolmogorov, Ito, and others (ref. 4). 

It might be supposed that a stochastic differential equation is any 
differential equation containing random terms or coefficients. However, the 
body of theory referred to above is more restrictive in its definition: 
randomness may only be introduced via a white-noise process. Our review 
starts naturally with a discussion of white noise and its integral, the 
Wiener process. 


The Wiener Process 

Consider the proverbial "drunkard's walk," in which the drunkard moves 
in one dimension by independent random steps, the n^*^ such step being 
denoted aW^^t* Each step is taken after a fixed time interval at. After 
N steps his position is 


N 

'^Nat " ^’^nat 
n=l 


In addition to being independent of one another, the drunkard's steps all 
have the same statistical properties. For modeling random dynamical systems 
one would like to replace this discrete random walk by a continuous random 
process. 

A continuous random walk might be constructed by introducing a continu- 
ous independent variable t = N at and letting N ~ while at dt + 0. 

If the resulting Wt is to be continuous, it must further be required that 
the steps become infinitesimal: aWp^t dW^, and of course, the above 

equation is now W^ = J' dW^,. With these stipulations W^ may be called 

the continuous drunkard's walk, although it is known more formally as the 
Wiener process. 

Since W-t is the sum of a large number of identically distributed 
independent random variables, by the central limit theorem (ref. 8) it has a 
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Gaussian probability density. The Gaussian distribution is fully character- 
ized by its mean and its variance . Our drunkard will not be allowed any 
preference in the direction of his walk, so has mean zero. 

Now, we return to the discrete drunkard's walk and compute its variance, 

~2 

'^Nat* steps of the drunkard's walk are statistically independent, so 

2 2 

^'^nAt ^'^mAt “ n ^ m and has a value, say a^, which is inde- 

pendent of n since the steps are statistically identical. (The overbar 
here denotes an average over an ensemble of realizations of the random walk, 
see equation (1-4).) Thus 



(II-l) 


9 

Remember that in the continuous walk, N At becomes t. But 

function only of t, and this will be the case in equation (II-l) 

? 

IS proportional to At; indeed, after suitable normalization a 

A 

equal to At and equation (II-l) tends to 


can be a 

, .. 2 
only if 

can be set 



(II-2) 


as N -*■ «> and At -► dt. Thus is Gaussian with mean zero and variance 

as given in equation (II-2), so its pdf is 

-w^/2t 

P(w.,t) = — (II-3) 

(ref. 8, p. 107; see the introduction for a description of the use of upper- 
and lower-case letters). 

The pdf P(w^) is required to be Gaussian by the central limit theorem, 

irrespective of the distribution of the increments dW^. However, it is 

natural to require P(dw^) also to be Gaussian; clearly, with mean zero. The 

variance of dW. is known too because dW. is the limit of aW ^ and the 
t L nAt 

2 

latter has variance a, = At. Thus (dW.) = dt and 

A t 


-(dw^)^/2dt 

P(dw. ) = -^ (II-4) 

^ V^27Tt 
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By their derivation equations (II-3) and (II-4) are consistent with 



As an aside, note that although is continuous in t, it strictly 

is not differentiable with respect to t. This can be seen from an order-of- 

magnitude analysis: dW^ = dt, so dW^ = 0(dt^^^) and dW^/dt = 0(dt“^^^). 

Therefore dW^/dt becomes infinite (almost surely) as dt -*■ 0. However, 

white noise dW-^/dt can be defined as a generalized stochastic process: it 

gets its name from the fact that its covariance dW^ dW^,/dt dt' is a delta 

function of t - t'. Although the incremental Wiener process dW^ divided 

by dt is usually called white noise, dW^ itself may be thought of as a 
white-noise process; of course, its variance is a factor of dt^ smaller 
than the usual white noise. 

Equation (II-2) gives the variance of the Wiener process at a particular 
time t. The two-time covariance is also readily determined from the proper- 
ties of W^. Because the increments in our drunkard's walk were assumed to 
be uncorrelated, W^ has independent increments; that is, 

fo t ^ t 

dW^ dW^ = < (1 1-5) 

dt T = t 


(Since dW-^ is Gaussian, equation (II-5) ensures complete independence of 
dW-{. and dW^. , x ^ t.) Covariances of W-^ follow directly from equation 
(II-5) and the definition of W^ as the integral of dW^: 


“A ■ I f ““f 

0 0 


min(t,T ) 


dW^. dW^„ 


= / 


min(t,T ) 


dt' = min(t,x ) 


(II-6) 


The first step uses the fact that ^ dW^, = W^ and dW^, = W^^^ - W^ 

are uncorrelated when a 0. This follows directly from equation (II-5); 
indeed, it is an alternative statement of the fact that the increments of 
which Wt is composed are independent. Multiple-time moments can similarly 
be determined from equations (II-4) and (II-5). 

The property of being composed of independent increments is a crucial 
property of the Wiener process. This makes it a continuous Markov process - 
our next topic. 


7 



Continuous Markov Processes 




A Markov process can be described as a time-dependent stochastic pro- 
cess in which the future is determined by the present, independently of the 
past . In deterministic dynamical systems this is the natural state of af- 
fairs; but a random system usually has a finite correlation time scale, which 
is the memory time of the system. Thus it is only in the special case where 
this correlation time becomes negligible that one obtains a Markov process. 

A continuous Markov process is also called a diffusion process (for reasons 
that should become clear later). Thus when describing turbulent transport 
we only refer to turbulent diffusion when we have in mind a Markovian model 
of fluid trajectories (eddy diffusion), but in general we refer to turbulent 
dispersion . 

The Wiener process is Markovian, for given its present value as W^, its 
future is determined as = W^ + dW^, where the integral is inde- 

pendent of the past. The objective of the present study of continuous sto- 
chastic processes is to model the random trajectories of fluid particles in 
turbulent flow. If we call the trajectory of a fluid particle Y(t), the 
trajectory generated by the Wiener process is the solution to 


dY/dt = dW^/dt 


(II-7) 


subject to an initial condition on Y. Clearly, the correlation time of the 
random velocity in equation (II-7) is zero (see eq. (II-5)), and this makes 
Y(t) a Markov process. Since dW-^/dt does not really exist, equations 
like (1 1-7) will hereinafter be written 


dY(t) = dW^ 


The random velocity on the right side of equation (II-7) is a function 
of time but not of particle position. Thus it describes diffusion in a homo- 
geneous medium. In an inhomogeneous medium a natural generalization of equa- 
tion (II-7) would be 


dY(t) = a[Y(t)]dW^ 


(II-8) 


where a(*) is some prescribed function of particle position and without 
loss of generality a(y) 0 for all y. Now the random velocity varies 
with position. (We will not consider unsteady media, where a(’) is also 
a function of t, although most of our discussion would remain valid.) 

Unfortunately, the interpretation of equation (II-8) is ambiguous 
unless the idea of a nonanticipating function is introduced. To see this, 
consider finite-difference approximations to equation (II-8). Compare 

Vl ' ''n " 


to 
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- ''n " - "n.t] 


(II-9b) 


observing that Y^ depends only on m £ n. Thus the right side 


of equation (II-9a) depends only on m ^ n, and on the increment 

W, - W But the Wiener process has independent increments, so 

(n+l)At nat 

'^{n+l)at ■ ^nat independent of W^^^, m <_ n. Therefore the statistical 
average of a(Y^) [W(r\+l)^t “ ensemble of Wiener processes is 


LW(„n)At - «nstJ - W L“(„n)At - ® ‘ ^ 

so equation (II-9a) gives rise to no mean particle veloc- 

ity. The second form does not have this property. In fact, if the truncated 
Taylor expansion 

- r -) ■ =C'n) * °<Vl - 

is substituted, equation (II-9b) can be solved for Y^+i : 

''nn = ''n * [“(n+D.t ' “nitl 


-I a(Y„)a'(V„) 


°(^n+l)at - ^nAt^' 


(II-IO) 


2 

Averaging equation (II-lO), using [W^^+2)At ~ '^nAt^ ^ gives 

\n - \ *7 ="'n) 

Thus equation (II-9b) gives rise to the mean velocity 

dY/dt = ^ a(Y)a'(Y) 

on letting At + dt . 
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In the continuous limit of equation (II-9a), a[Y(t)] is independent of 
dW-r, t ^ function defined in this way is called nonanticipating because 

it depends neither on future values of nor on the increment dW^. This 
is a convenient way to define a[Y(t)], and it will be adopted hereinafter. 
Practically, this has the effect that any function of Y(t) is independent 
of dWt, as in equation (II-9a). Thus equation (II-8) has zero mean 
velocity: 


dYlty = aLY(t)JdW^ = aLY(t)J x dW^ = 0 (II-ll) 

In general, particles in nonhomogeneous turbulence can have nonzero mean 
velocity. Such a velocity can be added to equation (II-8) to get the general 
form for a (one dimensional) diffusion process 

dY = b(Y)dt + a(Y)dW^ (11-12) 


By equation (II-ll), equation (11-12) has average 

dY = bTYldt 


Thus b( Y) is the mean particle velocity. 


Fokker-Planck Equation 

A diffusion process of the form of equation (11-12) describes the random 
trajectories of fluid particles through turbulent flow. The description of 
fluid motion in terms of particle trajectories is called a Lagrangian des- 
cription, as mentioned in chapter I. The alternative description discussed 
there is called Eulerian. The Eulerian analysis is based on a partial dif- 
ferential equation that governs the evolution of the probability density 
function (pdf) for the random process (eq. (11-12)). The pdf is a function 
of position in the fluid and of time and determines the probability of a 
trajectory passing infinitesimally close to a given point. Or, as an alter- 
native physical description, if a large number of marked particles are ini- 
tially released into the fluid, the pdf gives the fraction of these that at 
some later time will lie in a small volume around a given point. 

We now turn our attention to a derivation of the partial differential 
equation governing the pdf. In this equation y and t are the independent 
variables. Hence, it becomes necessary to define an average that can be 
evaluated at an arbitrary point in the fluid. Such an average is called a 
conditional average: we wish to calculate the moments of the increment 

dY(t) on the condition that only trajectories passing through y at time t 
are considered. However, once subject to Y(t) = y, the only randomness 
remaining in equation (11-12) is dW^; hence averaging is carried out 

readily. Making use of dW^ = 0, (dW^)^ = dt, and (dW^)*^ = O(dt)'^^^, we get 


10 


II III I I 


I III 


nil II I 


(dY)2 


dY = b(y)dt + a(y)dW^ = b(y)dt 


= (b(y)dt + a(y)dW^)^ = a^(y)dt + O(dt^) 


(dY)"^ = 0(dt^^^) , n > 2 




(11-13) 


For later purposes it is only necessary to keep the terms in equation (11-13) 
that are of order dt. 

By definition, the future of a Markov process is determined by its 
present state alone. Thus the probability density at time t is determined 
by that at time t - dt and by a transition probability: 

P(y,t) = / P(y - dY, t - dt)dP.(dY; y - dY) (11-14) 

all ' 

dY 


Here P(*,*) is the pdf for fluid particle positions and dPj(*,*) is the 
transition probability that a particle moves from y - dY to y in time 
dt. The integral is over all possible jumps dY. The Fokker-Pl anck equation 
is derived from equation (11-14) as follows: The right side of equation 
(11-14) is expanded in a Taylor series in dY and dt, and the following 
definitions of the moments of dY: 


f dPy (dY;y) 





dPy(dY;y) 


(dY)"(y) 


\ (11-15) 



8 dP-^(dY;y) 

sy 



dPj(dY,y) = 


3(dY)^y) 

3y 


are then substituted. The (dY)>^ are given explicitly by equation (11-13). 
If they are substituted, keeping only terms up to order dt, one finds from 
equation (11-14) 


P(y,t) = P(y,t) 


_ 9P(y _ , _ t) _ + 1 3^P(y . »t„) , 

3 ^ 2 3^2 


dt 


3P(y»t) 

3t 


3 dY 

ay 


P(y,t) 


= P(y.t) + dt 


1 3^[a(y)P(y,t)] _ 3P(.y,t) _ 3[b(y)P(y,t)] 

2 ...2 3t 3y 


+ O(dt^) 
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The last equation can be written as 


3P(y,t) ^ 3[b(y)P(y,t)1 1 3^[a(y)P(y,t)] , 

3t ay - 2 3^2 ^ 

Equation (11-16) is called the Fokker-Planck equation for equation 
(11-12); clearly, it is a diffusion equation for P(y,t). Indeed, if 
3^(y)/2 is defined as the eddy diffusivity, a^(y)/2 = K(y), then equa- 
tion (11-16) can be put into the form 

— + — T/^b - — ^ P 

3t ay ayy 

which is the conservation equation form of an eddy diffusion equation. In 
equation (11-17) [b(y) - aK(y)/ay] plays the role of an Eulerian mean con- 
vection velocity. When the latter is zero, b(y) = aK(y)/ay and equation 
(11-17) becomes 

iZ = i_ /k — ^ 

at ay \ ay/ 

with the corresponding form of equation (11-12) being 


dY = (Y)dt + y/2 K(Y)dW^ (11-18) 


-hH) 


(11-17) 


(since a^ = 2K and b = aK/ay). Averaging equation (11-18) over an 
ensemble of trajectories gives the Lagrangian mean velocity 


^ _ IFm 

dt ay ' ' 


(11-19) 


Hence, in a nonhomogeneous medium (where K is a function of y) the 
Lagrangian mean velocity of fluid particles may be nonzero, even though their 
Eulerian mean velocity is zero. More will be said about equations (11-18) 
and (11-19) in chapter III. 


Ito's Theorem 

We return from our excursion into Eulerian analysis to equation (11-12) 
and investigate its calculus. Because the functions in equation (11-12) were 
defined to be nonanticipating, certain peculiarities arise; it is necessary 
to derive a calculus consistent with the nonanticipating property. 

Consider the following example: In equation (11-12) let b = 0 and 

a = Y, 

dY = Y dW^ (11-20) 


One is inclined to integrate this as 


12 



I 


Y = Ce *' (11-21) 

where C = constant. However, on averaging equation (11-20), using the non- 
anticipating property to write Y dW^ = ^ = Y x 0, one finds dT = 0 or 
T = C. But by equations (II-3) and (1-5), equation (11-21) has the average 




Ce 


t/2 


Hence equation (11-21) is not a nonanticipating solution to equation (11-20) 
Clearly Y = C exp(Wt - t/2) has the correct mean, and in fact it is the non- 
anticipating solution to equation (11-20). How does one know this? 

The derivative of a function f(Wt) of Wt can be defined as df 
in the relation 

f(W^,t) = J' df(W^,t) + constant 


Expanding the integrand in a Taylor series, to order dt, gives 

A-/ dt*o(dt) 


It can be shown that (ref. 3, p. 98) 




dt 


(11-22) 


(11-23) 


using equality in a statistical sense. (Roughly, the equality (11-23) is 
shown as follows: In mean square, upon averaging over dW^ with fixed 



f 

JJ < 

■If 


( dW ^)2 ( dW ^,)2 


dt (dW^,)‘ 


3^f g2x 

(dt df - 2dt 

3W^ 3W^, 


dt df 

df + dt df ) = 0 

(11-24) 
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t t' . 


Since the 


having used (dW^)^ = dt and (dW^) (dW^,)^ = dt dt‘, 

square of a number is nonnegative, the fact that the mean value of the left 
side of equation (11-24) is zero implies equation (11-23) (except on a set 
of measure zero. ) 

Substituting equation (11-23) into equation (11-22) and equating quan- 
tities under the integrals in that equation give 


8f(W ,t) 


j 3^f(W^,t) 3f(W^,t) 


3W, 


3t 


dt 


(11-25) 


As an example of the formula (11-25) for differentiation, if 
f = exp(W^ - t/2) 


df = f dW^ + -^ f dt - ^ f dt = f 


dW, 


This shows exp(W^ - t/2) to be a solution to equation (11-20), as was 
claimed previously. 

The formula analogous to equation (11-25) for functions of Y(t) is 
called Ito's theorem. If in equation (11-22) the function is f[Y(t),t], 
then 




o(dt) 


/ 


If •><'')<'* ^ W ^ a2(Y)(dW^)2 + If dt 

3 I 


+ o(dt) 


having used equation (11-12) for dY. Again, we use equation (11-23) and 
equate integrands to find 


df = 


b(Y) 


3f 

3Y 


.1 


a^(V) 


3^f 

3Y^ 


3f 

3t 


dt + a(Y) 


3f 

3Y 


dW. 


(11-26) 


This is Ito's theorem. 


Attainability of Boundaries 

As an application of theorem (11-26), consider a turbulent flow between 
boundaries at y+ and y_. Consider a trajectory that originates at some 
point yQ in the fluid and strikes one of the boundaries at some arbitrary 
later time. We will determine the probability that this trajectory strikes 
the boundary y+ (or y_) . If this probability is zero for one of the bounda- 
ries, that boundary (almost surely) cannot be reached from within the fluid 
and is called nonattainable. 
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Let 



(11-27) 


where f(y) is the solution to d^f/dy^ + (2b/a2) df/dy = 0, satisfying 
f(yo) = 0; tl^us by equation (11-26) 


df(Y) = a(Y) exp 



b(x) 

a^(x) 


dx 


dWt 


Averaging this gives df [Y(t) ] = 0 (cf. eq. (II-ll)); hence f[Y(t)] is a 

constant. But f(yQ) = 0, so f [Y(t) ] = 0 independently of t. 

Consider trajectories that at some time strike a boundary y+ or y_. 
There is a probability P+ that a given trajectory strikes y+, and P_ 
that it strikes y_. These are the probabilities that we wish to determine. 

Since f ( Y) =0 is true at any time, it is true at the particular time when 
the boundary is struck. At that time, Y is at either y+ or y_, so 

fW = f(y+)P+ + f(y_)P_ = 0- Or, because P+ + P_ = 1, 


- f(y+) 

= Lf(y±) - f(y+)] 


(11-28) 


Equation (11-28) gives P± explicitly in terms of the integral (11-27) with 
a nonrandom upper limit of y+ or y_. In particular if the integral does not 
converge at y_ (say), then f(y_) = and P_ = 0 and thus y_ is non- 

attainable. 

In the case of equation (11-18), a^(y) = 2K(y) and b(y) = K'(y), and 
equation (11-27) becomes simply 


f(y) 



K(yo) 

Kirr 


dy ' 


Then f(y_) is infinite if K tends linearly to zero at y_. Such a 
boundary has P_ = 0 and it cannot be struck by a typical fluid trajec- 
tory. Generally, if b(y) K'(y) and K(y) +0 as y -► y_, y_ is 

not attainable. In the inhomogeneous turbulent flow near a surface it is 
usually the case that K 0 at the surface. 
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Numerical Methods for Stochastic Differential Equations 

Often closed-form solutions to equation (11-12) cannot be found and one 
wants to solve it numerically. This requires replacing dt by at and dW^ 
by a Gaussian random variable with mean zero and variance at. Two consid- 
erations present themselves in connection with approximating equation (11-12) 

1/2 

consistently: first, dW. becomes 0(at ' ) so that an approximation to 

0(at”) requires retaining terms to 0(dW^^); second, the approximation must 

preserve the nonanticipating property of the coefficients a(Y) and b(Y). 

Consistent finite-difference approximations to Ito differential equa- 
tions are derived in references 9 and 10; here the method of reference 9 is 
followed. This method consists of integrating equation (11-12) between times 
t,^ and tp+i, separated by at 

Vl = ^n [b(Y)dt + a(Y)dW^] (11-29) 

and approximating the integral by successive approximations. Thus, to 
0(atl/2) 


'n+1 


= Y. 


a(Y^)aW 


where aW = is a (computer generated) Gaussian random vari- 

2 

able with aW = 0 and (aW) = at. The preceding first approximation can be 
generalized to Y(t) = + a(Y^)(W^ - W^^^) ; t > t^. Using this in equa- 

tion (11-29), to 0(at), yields 

Vl - ''n " ^ 


. Y„ - a(Y„)AW * b(Y„)at * a(Y„)a'(Y„) («, - dW^ 

= Y„ * a(Y„)aW * b(Y„)at * a(Y„)a'(Y„) W^dW^ - 


(11-30) 


Now, according to Ito's theorem (eq. (11-26)) 


dw; .. 

= W, dW, - ^ 
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and thus 



"t "“t = I ("(nn)At - “Lt - 



Substituting this into equation (11-30) gives, to 0(at), 
Ynn = Y, - a(Y^)AW - b(Yjat - 


a(Y )a'(Y ) . 

— [(aW)*^ - At] (11-31) 


Mil'shtein's (ref. 10) method for deriving equation (11-31) consists of 
showing that it approximates Yp+i with a mean square error of 0(At3). 
Higher order finite difference schemes are derived in references 9 and 10. 
Mil'shtein also notes that it may often be convenient to replace equation 
(11-31) by a Runge-Kutta scheme with the same order of accuracy: 

''nn - ''n * " ^1 

+ P 2 spn * “1 * “2 (n-32) 

where pi + P 2 = 1 and P 2 c»i = -P2«2 = 1/2. Expanding the last term to 
0(At) recovers equation (11-31). Formula (11-32) has the advantage of not 
requiring a derivative of a(Y). 


Types of Boundary Conditions 

When dispersion takes place in a bounded fluid, one must impose boundary 
conditions, and these conditions must be compatible with the stochastic dif- 
ferential equation. Numerical formulas such as (11-31) and (11-32) also must 
be supplemented with boundary conditions. 

Reference 4 classifies the types of boundary that can occur and the 
boundary conditions compatible with them; the appendix to this chapter des- 
cribes this classification. In the present context the most important con- 
sideration is whether or not the boundary is attainable (eqs. (11-27) and 
(11-28)). If a boundary is not attainable, no boundary condition can be im- 
posed on it; stated otherwise, it automatically satisfies a no-normal-flux 
condition. An attainable boundary will usually be refl ec ting or absorbing . 
Other variants could be conceived: for instance a trajectory could remain 
at the boundary some time before leaving, modeling an initially absorbing 
boundary that saturates; or the boundary could be an interface to a region 
with different properties, modeling an impedance condition. Here only re- 
flecting, absorbing, and nonattainable boundaries are considered. 


Nonattainable Boundaries 

Although in principle these boundaries cannot be reached, when a finite- 
difference approximation is made, it becomes possible for a trajectory to 
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reach a nonattainable boundary. Such a trajectory must return to the in- 
terior in a suitable way. The following scheme might be used in a numerical 
calculation to handle this eventuality. 

If on one time step a trajectory reaches the "nonattainable" boundary 
y = 0, at the next time step it may be considered a trajectory originating 
at y = 0. Its position then can be chosen from the probability distribution 
for particles originating on the boundary; the relevant pdf is derived next. 
(To avoid confusion, it should be remarked that a nonattainable boundary can 
be the origin of particle trajectories. Its property is that once these 
trajectories move into the interior, they cannot return to the boundary.) 

Suppose K(y) -*■ ay and b(y) + 3 near the boundary, where a 
and 3 are constants with 3 2 “ 6 < a, the boundary is attainable). 

Then equation (11-17) becomes 


^4. 

at 




(yP) 


The pdf at time At for trajectories originating at the boundary y = 0 is 
the solution to this at t = At, subject to P(y,t = 0) = 6(y): 


P(y,At) 


( 3-a) / o 



-y/aAt 

6 

r ( 3 /a ) a At 


(11-33) 


where it is assumed that a At << 1. (r(*) here is the gamma, or fac- 

torial , function. ) 

In a numerical calculation, if a trajectory reaches the boundary on some 
time step, its position at the next time step should be a random variable 
Y chosen from the distribution (11-33). At subsequent time steps its posi- 
tion is determined by equation (11-31) or by equation (11-32), until the 
boundary is again reached. 

Next, we consider attainable boundaries at which K(y) ^ 0. This is the 
case where the turbulent eddy diffusivity is nonzero at the boundary. 


Reflecting Boundaries 

If the boundary is at y = 0 and the flow domain is in the region 
y > 0, reflection is achieved by replacing Y by IyI as soon as Y be- 
comes negative. Equivalently, the condition of reflection is imposed by 
placing a particle, which at the end of a time step has passed a distance 
through the boundary, an equal distance above the boundary. 

An analysis of the reflection boundary condition proceeds through con- 
sideration of the case where the flow domain is the entire region y > 0. 

In this case a solution Y(t) to 

dY = b(Y)dt + a(Y)dW^ (11-34) 

satisfying Y(t) > 0, and being reflected at y = 0 is required. It will be 
assumed that a(0) > 0. Note that as formulated, equation (11-34) governs 
Y(t) only when Y(t) > 0; Y(t) at y = 0 is governed by the boundary con- 
dition. What we will do is find another process Y'*’(t) defined on the entire 
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y-axis, and hence not subject to any boundary condition, but which satisfies 
equation (11-34) in y > 0. The process T{t) also will have the property 
that lY^(t)l satisfies equation (11-34) for all nonzero values of Y'^(t) ; 
|Y'*’| satisfies the reflection condition and hence is the required solution 
Y(t) = |Y+(t)| . 

In equation (11-34) b(Y) and a(Y) are defined only for Y > 0. Their 
definitions will be extended to all Y by letting b(-Y) = -b(Y) and a(-Y) 
= a(Y); that is, b(Y) = sgn(Y)b(|Yl) and a(Y) = a(|Y|). It will be seen 
that these extended definitions give Y'*’ the properties required of it in 
the last paragraph. 

Let Y‘‘‘(t) satisfy 


dY"^ = b(Y‘*’) + a(Y‘^)dwJ 


(II-35a) 


Then by Ito's theorem (eq. (11-26)) 


dlY^I = b(Y^) sgn(Y^)dt + a(Y^)sgn(Y^)dW^ + ^ fi(Y^)a^(0)dt 


But dW. and -dW. are statistically equivalent, so we may replace 
sgn(Y )dW^ by dW^ and use the extended definitions of a(Y) and b(Y) to 
rewrite this last equation as 

dlY'"l = b( |Y'"l)dt + adY'^DdW^ + | 6( 1 Y^| )a^(0)dt (II-35b) 

Thus 1 Y^l satisfies all of the criteria that were set for it: it is deter- 

mined by the equation (II-35a), on which no boundary conditions need be im- 
posed, and it is a solution to equation (11-34) irrespective of whether 
Y'‘'(t) is positive or negative. Thus, Y(t) = |Y‘'’(t)l. 

Equation (II-35b) shows that when the boundary condition is incorporated 
into the dynamical equation, it takes the form of an impulsive velocity at 
the boundary - the 6-function term. This impulsive velocity accomplishes 
the reflection of the trajectory. It will be shown below equation (11-43) 
that a consequence of the 6-function in equation (II-35b) is to make the pdf 
for Y satisfy the boundary conditions 3P(y,t)/3y =0 on y = 0. 

The reflection condition is illustrated clearly by the simple case 
b = 0, a = o = constant. In this case equation (II-35a) is easily inte- 
grated to 

y" = yg - aw^ 


where yg is the integration constant. Therefore Y = lY^I = lyg + A 

given value of Y occurs when yg + aW^ = Y and when yg + aW^ = -Y. 
Therefore the probability (density) of observing a given Y is 


P(y,t) 


= P 


(«t 




(11-36) 
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-wf/2t 


But P(w^) = e ^ lyfHt (eq. (II-3)), so 


r 


P(y,t) = 


a •y^2iTt 


< exp 


(y - yg)' 

2a 


+ exp 




(y + yg) 


2a"t 


^ >H(y) (11-37) 


-IJ 


(The Heaviside function H(y) times 1/a arises on transforming from 
wt to y = ly"^! as the independent variable.) Equation (11-37) is the pdf 
for y, subjected to a reflection boundary condition, given that, at t = 0, 

y = yo- 


Absorbing Boundaries 

It is shown in the appendix to this chapter that some boundaries are 
naturally absorbing, in the sense that once a trajectory reaches such a 
boundary it cannot escape from it. If the boundary is not naturally absorb- 
ing but is attainable, one can still impose the condition that, once the 
boundary is reached, the trajectory stops. As might be expected, the formal 
analysis of absorption (provided by ref. 4) shows that the absorption condi- 
tion corresponds to imposing the boundary condition P(y,t) =0 at y = 0 
on the Fokker-Planck equation (11-16). This is discussed further in the next 
section. 

Again this boundary condition is readily illustrated through the case 
b = 0, a = a. Since trajectories terminate upon hitting the boundary, the 
probability (density) of observing a particle at y at time t is the 
probability of observing a particle at y in the absence of a boundary, 
minus the probability of the particle's trajectory reaching y = 0 at some 
time prior to t. Now, dW^ is statistically equivalent to -dW^, and dW^ 
is uncorrelated with dW^ 4 r if 4 o. Thus if we take a trajectory 
dY = o dW-f- that reaches y at time t, having crossed y = 0 at some prior 
time, and from the first time that trajectory reaches y = 0 onward replace 
dW|- by -dW^, we obtain an equally likely trajectory that arrives at -y at 
time t. Hence there is a one-to-one correspondence between trajectories that 
reach y having crossed y = 0 and trajectories that reach -y. This 
correspondence can be used to derive the pdf for Y(t). Thus as described 

P[Y(t) = y I absorbing boundary] 

= P[Y(t) = y 1 no boundary] - P[Y(t) = y | no boundary and 
Y(t') = 0 for some t' < t] 

= P[Y(t) = y I no boundary] - P[Y(t) = -y | no boundary] 

(11-38) 

Here the vertical bars are to be read "given the condition to the right." 

Since Y(t) - yg = aW^ and P(W-t,t) is given by equation (II-3), 

P[Y(t) = y I no boundary] is 

exp[^-(y - yQ)^/2a^t]/a -y/T^rt 
With this, equation (11-38) is 
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P(y,t) = 


1 


yQ)^/2a^t] - exp[-(y + yQ)^/2a^t])H(y) (11-39) 



In equation (11-39) P(y = 0,t) =0 as expected from the first paragraph of 
this section. 


Application to Mean Concentration 

Having now summarized the theory of continuous Markov processes, it is 
time to return to our starting point, equation (I-l), and see how all of this 
relates to it. The mean concentration of contaminant at y,t is determined 
from any initial distribution of concentration by equation (I-l). The aver- 
aging in this equation is over random initial positions Y(0), given the 
final position Y(t) = y. In the overbar notation of equation (1-4), equa- 
tion (I-l) is 


C(y,t) = CqLy(0)J 


(11-40) 


Thus calculation of C(y,t) requires trajectories to be integrated backward 

in time from their known position y at time t to their random position 
Y(0) at time zero. The known initial concentration is then evaluated at the 
random initial position, and C(y,t) is the average of all such evaluations, 
as in equation (11-40). 

To calculate the requisite reversed trajectories, the reversed time s 
is introduced: ds = -dt. As time runs from 0 to t, s runs from t to 

0. In terms of s, equation (11-12) becomes 


dY(s) = -b[Y(s)]ds + a[Y(s)]dW^ 


(11-41) 


with Y(0) = y. We next show that the calculation involving this reversed 
diffusion does indeed give the solution to the actual forward diffusion 
process . 

Let Co(y,s) satisfy 


aCQ(y,s) 

Ts 


b(y) 


3Co(y,s) ^ 1 2, V 

ay 2 ® 


(y,s) 


(11-42) 


with initial condition CQ(y,s = 0) =T(y,t). (Recall that s = 0 corresponds 
to time t.) Then by Ito's theorem, substituting -3/3s for 3/3t, we get 


aCo(Y,s) 

dCQ(Y,s) = a(Y) dW^ 

Integrating this between s = 0 and t yields 
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a(Y) 


(11-43) 


CQ(y,s = 0) = CQ(y,s 



sCqCY.s) 

iY 


dW 


s 


After averaging and reverting to forward time, so that s = 0 is replaced 
by time t and s = t is replaced by time 0, equation (11-43) becomes 

r(y,t) = Cq[Y( 0)] , which is equation (11-40). Here the initial condition 
CQ(y,s = 0) = C(y,t) is used. 

Thus equation (11-40) is obtained provided that Cq satisfies equa- 
tion (11-42). Equation (11-42) is quite similar to the Fokker-Planck equa- 
tion (11-16). Indeed, it is called the backward Fokker-Planck equation 
(ref. 3, p. 159). (Note that we have let s increase in the direction of 
negative t; usually ds is taken as positive in the direction of forward 
time so the first term in eq. (11-42) is -aCo/as.) It is easily shown that 
this equation with fixed final time t and variable initial time s corre- 
sponds to exactly the same diffusion process as does equation (11-16), where 
the initial time is fixed and the final time is variable (refs. 3 and 4). 
Therefore the correctness of equation (11-40) as calculated by reversed dif- 
fusion has been demonstrated. 

As an addendum, it is now clear why, when a reflecting boundary is pres- 
ent, Co(y»s) must satisfy aCo(y,s)/3y =0 on the boundary. The delta func- 
tion in equation (II-35b) adds a term 


1 



6(Y)dt 


to the right side of equation (11-43). In deriving equation (II-35b), it was 
assumed that a(0) 4 0, so aCQ(0,s)/3y must equal zero if equation (11-40) 
is to obtain. Thus, calculation of T(y,t) from the Fokker-Planck equation 
with a zero normal derivative gives the same result as expression (11-40) for 
reflected trajectories. 

A similar line of reasoning can be followed in the case of absorbing 
boundaries. Let Co(y,s) be as above and satisfy equation (11-42). Let 
yU) be the solution to equation (11-41), now subjected to the constraint 
that if, at any s < t, Y reaches the boundary y = 0, the trajectory sub- 
sequently remains there. For trajectories that do not reach the boundary, 
equation (11-43) will apply; but for a trajectory that reaches the boundary 
and remains there, Y(s = t) = 0, with the effect that equation (11-43) 
becomes 

/ 3 C.(Y,s) 
a(Y) — 

u 
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Consequently, if the indicator function 


X 


if the trajectory reaches the boundary 


otherwise 


is introduced, averaging equation (11-43) gives 


C{y,t) = (1 - x)Cq(0,s = t) + xCq[Y(s = t),s = t] (11-44) 

Only if Cq(0) = 0, so that the first term in equation (11-44) vanishes, will 
equation (11-40) obtain. (It should be realized that the factor of x mul- 
tiplying the second term in eq. (11-44) only indicates that averaging is 
carried out over those trajectories that are not absorbed at the boundary.) 
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APPENDIX: CLASSIFICATION OF BOUNDARIES 


Gikhman and Skorohod (ref. 4, chapter 5) classify four types of bound- 
ary, according to the convergence of 




where the boundary i s at y = y_, trajectories originate at y = yo > y_ 
and y_^ > y^ is an arbitrary interior point of the fluid. Their classifi- 
cations of the boundary are 

(1) Nonattainable: = “ 

(2) Asymptotic: < “, ^2 ~ 

(3) Absorbing: < ”, 1-2 = “ 

(4) Normal: < “, L 2 < ”, < “ 

A nonattainable boundary is one that cannot be reached from the interior 
of the fluid. An asymptotic boundary is one that trajectories can only reach 
asymptotically as t 0 °. An absorbing boundary is one that is attainable 
but from which trajectories cannot escape in finite time. A normal boundary 
is attainable and not absorbing. Clearly, no boundary conditions can be im- 
posed on boundaries of types 1 or 2. Continuous trajectories can only sat- 
isfy an absorption boundary condition on type 3 although, since this boundary 
can be reached in finite time, one can impose the condition that, having 
reached the boundary, trajectories jump from it to a random position within 
the fluid. On boundaries of type 4 one can arbitrarily impose a condition 
such as reflection or absorption. We will give a brief derivation of this 
boundary classification. 
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The function is equivalent to the function f(y+) - f(y_) in equa- 

tion (11-28); the fact that = * implies that the boundary is unattain- 
able is explained in the main text. 

To explain type 2 and 3 boundaries, it is necessary to introduce the 
idea of a stopping time t . This is simply the (random) time at which a 

a trajectory originating at^ yo reaches a boundary of the interval (y_,y+). 
Applying Ito's theorem (eq. (11-26)) to 


J' ^ 9 ( X ) dx 

f(y) - 


/•y+ 

/ cp(x)dx 


/?' c 


dz 


<p(z)a (z) 


dx 


- 2 


/ b t 


dz 


cp(z)a (z) 


dx (II-A3) 


where 


tp(x) = exp - 


vf 


b(z) 

a^(z) 


dz 


gives 


df = -dt + a(Y) 


df(Y) 

dY 


dWt 


Integrating from 
and that f(y+) = 


t = 0 to t = T 

•^0 

f(Y_) = 0, gives 


noting that 



by definition 


f(yo) 



df(Y) 

dY 


dW, 


(II-A4) 


Averaging this gives 


f(yo)-fy^ (n-A6) 

Since y+ is an arbitrary point, it generally can be reached in a finite 
time. Hence can be infinite only if (on average) it takes an infinite 

time for trajectories to reach y = y_. (Ref. 4 shows that "on average" in 
the last sentence can be changed to "with probability 1.") Thus f(yo) = “ 
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implies that the boundary is of type 2. But it is clear from equation 
(II-A3) that f(yo) is infinite if and only if L 2 is infinite. ^ 

If a boundary is attainable in finite time, Li < “ and L 2 < “, a third 
possibility exists: having attained the boundary, a particle may be unable 

to escape it in finite time. To derive the condition for this to occur, we 
consider the stochastic process in which a particle is released at yg and 
each time it reaches the boundary of the interval (y_,y+) it immediately re- 
turns to yg. Then we let yg y_ and determine whether there is nonzero 
probability that after a finite time the particle will have reached y+. If 
this probability is zero, the particle (almost surely) cannot escape from 
y_ because y+ can be taken arbitrarily close to y_. 

After a sufficiently long time t, it is expected that the number of 
times the boundary has been reached in the preceding stochastic process will 
be 


N = t/T 


^0 


(II-A6) 


1 + 

The mean stopping time can be decomposed into the mean stopping times t for par- 

— -^0 ^0 

tides reaching y+ and t for those reaching y . (This does not seem to be discussed in 

" ^0 

ref* 4.) Let P+ be as in eq. (11-28). Then, the stopping time is decomposed as 


T = T P + T P 


V + 


^0 - 


The mean stopping time is computed as follows: Let 

^0 


g(y) = t <P(x)dx j ^ tp(x)dx 


y ' y 

with 9(x) defined below eq. (II-A3). Note that g(y+) = 0, g(y.) =1, and g(yo) = P_* Let 


h(y) = 2g(y) I 9(x) 




g(z) 


a"(z)<P(z) 


r 

dz|dx - 2 I <P(x) / 


g(z) 


dz 


a (z)cp(z) 


dx 


(For computing replace g(z) by 1 - g(z) in the innermost integral.) Note that 

h(y. ) = h(y ) = 0. Then by Ito's theorem 


d[h(Y) + tg(Y)] = a(Y)[h'(Y) + tg'(Y)]dW^ 


Integrating this from time zero to the stopping time and averaging give 


h(yo) = 



or 


^y ” * Since g(y ) = 1, the condition for convergence of ^i(yQ) is the same as that 

for eq. (II-A3), confirming that L^ = X means y can only be approached asymptotically. 
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where it is understood that N is the integer part of the right side of 
equation (II-A 6 ). The probability that the boundary struck is y+ on any 
one of these times is P. as given in equation (11-28). In N releases the 

N 

probability that y_ is struck every time is P_. Hence the probability 
that y^ is struck at least once is 

1 - p|^ = 1 - (1 - PJ*^ (II-A7) 


Now, by equations (II-A5) and (II-A3), as y^, tends to y , t and P^ 

tend to zero. Therefore N ~ and P^ 0 in equation (II-A7) and it 
becomes 


1 


- exp 


- t 



(II-A 8 ) 


This will be greater than zero only if 


1 im 





< “ 


A direct calculation using equation (II-A5) for 
for P+ yields 


t and equation (11-28) 
^0 



(II-A9) 


Since and were assumed to be finite, this last equation shows that 

equation (II-A 8 ) is nonzero if and only if L 3 < “. 
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CHAPTER III: EXAMPLES OF EDDY DIFFUSION 


In this chapter, three examples illustrating the statistical approach 
to turbulent dispersion theory are considered. These examples serve also to 
make the theory of the previous chapter more concrete. The examples consid- 
ered have been dealt with in various publications; however, they are analyzed 
here entirely through consideration of the solutions to stochastic ordinary 
differential equations, rather than the more usual treatments based on the 
diffusion equation. Because the Fokker-Planck equation relates stochastic 
ordinary differential equations to the diffusion equation, there clearly must 
be an equivalence between the present methods and those based on the diffu- 
sion equation. For the most part the second and third examples are treated 
more easily by analysis of the Fokker-Planck equation than by the approach 
used in this chapter. Our purpose here, however, is to show that statistical 
reasoning can be used equally in these examples. In doing so we hope to 
make clear that the often-made distinction between "statistical" theory and 
"eddy diffusion" theory is fictitious: eddy diffusion theory is really just 

a Markovian statistical theory; the eddy diffusion equation is just the 
Fokker-Planck equation for the pdf of a Markovian random process. The 
following examples show, furthermore, that by analyzing the random process 
itself one can obtain the same results obtainable from its Fokker-Planck 
equation. In fact, for example 1 the present analysis seems simpler than 
the diffusion equation analysis. 

Aside from their pedantic value the statistical methods used in this 
chapter will be useful in chapter IV, where a slightly more complex stochas- 
tic model of turbulent dispersion is considered. It is hoped these methods 
will also provide an insight that will be of use to others who wish to form- 
ulate and analyze stochastic models. 

As a further word of preface to these examples of eddy diffusion, it 
should be observed that the eddy diffusion approximation is only justified 
in an asymptotic sense. It applies at times that are long compared with the 
correlation time (the Lagrangian time scale) of the turbulence. This is 
discussed further in connection with example 2 and in chapter IV. 

In the following examples dispersion takes place in the y-direction and 
is governed by a stochastic equation of the form 

dY(t) = K'[Y(t)] + ^2K[Y(t)JdW^ (III-l) 


where K'(y) = dK(y)/dy. However, we will consider two-dimensional flows, 
with a mean shear in the x-direction. The turbulent velocity in this direc- 
tion will be ignored with respect to this mean velocity; thus 

dX(t) = U[Y(t)]dt (III-2) 

where U(y) is a known mean velocity, such as the log profile 
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(III-3) 


^ 1" (yj) 

In this, u* is the friction velocity and yg is the surface roughness. 

The Fokker-Planck equation for equations (III-l) and (III-2) is, in analogy 
to equation (11-16) , 

9P(x,y,t) ^ 3l!U(y)P(x,y,t) ] _ ^3 K(y) ^P(x,y,t) (III— ^ 

3t 3x 3y 3y 

(refs. 3 and 4). Equation (III-4), of course, is the diffusion equation 
usually used in analyses of the following examples. 


Example 1: Shear Dispersion 

In this first example we consider diffusion in homogeneous turbulence, 
so K in equation (III-l) is a constant: 

dY(t) = y/W dW^ (III-5) 

If the turbulence is homogeneous, to be consistent, any mean shear ought also 
to be homogeneous; that is, the rate of shear should be a constant. If this 
rate of shear is n, then U(y) = ny and equation (III-2) becomes 

dX(t) = nY(t)dt (III-6) 


This example of homogeneous turbulence in uniform shear might approxi- 
mately describe horizontal dispersion in a horizontally infinite boundary 
layer. For instance, it might represent dispersion in the atmosphere when 
horizontal variations of the wind are present. 

It is sufficiently general to use the initial condition Y(0) = 0. Then 
equation (III-5) integrates to 


Y(t) = (III-7) 

From the properties of (eqs. (II-2) and (II-3)) Y(t) is a Gaussian 

—7 ~2 

random variable with mean zero and variance Y = 2KW^= 2Kt: 

P(y,t) l^g-y^/4Kt (III-8) 

2-^irKt 


The function 
(III-6) 


X(t), subject to X(0) = 0, is found by integrating equation 



(III-9) 
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Since the sum of a number of Gaussian random variables is a Gaussian random 
variable, is Gaussian. Hence X is also Gaussian, and from the 
properties of (eq. (II-6)) 


X = 0 

X^ = 2n^K W. ,W. „ dt' dt" 

*'o *'o ^ ^ 


= 2n^K min (t',t") dt' 

•'o •'o 


2n^Kt^ 

3 





(III-IO) 


Thus 


P(x,t) = 


-3x^/4Kn^t^ 


^irKn^t^/3 


(III-ll) 


Clearly, X r^nd Y are jointly Gaussian random variables, so their 
joint distribution is given by reference 8 (eq. (5.1.28)). This joint 

distribution requires that we evaluate X(t)Y(t) : 


XY = 2KtiW. f W. , dt' 
t t 



dt' 


Knt^ 


(III-12) 


Substituting equations (III-IO) and (III-12) and Y = 2Kt into equation 
(5.1.28) of reference 8 gives 


P(x,y,t) 


2TrKnt^ 



3x^ + ^ _ 3xy 

Kn^t^ Knt^ 


(III-13) 


Physically, P(x,y,t) is the concentration distribution that would evolve 
from an instantaneous point release of contaminant at x = 0, y = 0, t = 0. 
Conditional x-moments . -i Expression (III-IO) is the second moment of 

this contaminant cloud in thej x-direction (cf. eq. (1-4)) 

1 


X^(t) = 



P(x,y,t)dx dy 


31 


It is a moment of the entire doud. Sometimes one is interested in the 
x-moments at fixed y, rather than in the moments of the entire cloud. These 
are usually calculated in terms of the conditional pdf, P(x|y): 


x^(y.t) 



P(x y,t)dx 


The conditional pdf is the pdf of observing a particle at a random point x, 
given that its other coordinate is y. This pdf is given by 

P(x|y,t) = (III-14) 


If written P(x,y) = P(x|y)P(y), equation (III-14) is intuitively clear: the 

unconditional pdf is just the probability of observing a particle at x, 
given its y-position, times the probability of observing a particle at that 
y-position. 

From equations (III-8), (III-13), and (III-14) 


P(x|y,t) = 




exp 


— ^ (x - nty/2)^ 


(III-15) 


That is, X here is a Gaussian random variable with mean nyt/2 and variance 
Kn^t3/6; these are the first and second moments of x at fixed y. 

Wiener process with fixed end points . - Another way the problem of de- 
termining moments of x at fixed y can be approached is through consider- 
ation of a Wiener process with given end points. That is, we consider the 
subset of all random trajectories W^,, 0 < t' < t that satisfy W« = 0 and 

t t 

W^ = y. This stochastic process will be denoted by W^,; thus Wq = 0 and 

W^ = y. The processes W^ and W^, are illustrated in figure III-l by 
schematic drawings of ensembles of trajectories. Because of the stationarity 



Figure III-l. - Schematic comparison of Wiener process and process with 
fixed end points W}i. In the former the cloud of trajectories spreads mono- 
tonically from the origin. In the latter the cloud spreads and then collapses 
about the mean trajectory. 
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of the increments of W^, on average W^, ought to move linearly between 0 

and y; that is, W^, = t'y/t. A stochastic process with this average, and 
which satisfies the end conditions is 


Because W^, and are Gaussian, this is a Gaussian random process with 

variance 



(III-17) 


Another consequence of the increments of being stationary is that 

t ^ 

the variance of W^, should be symmetric about t/2; that is, if t' is 

replaced by t - t', equation (III-17) retains its value. (To see that this 
must be so, consider the case where y = 0. Then trajectories originating at 
(y = 0, t* = 0) and ending at (y = 0, t' = t) must be statistically equiva- 
lent to trajectories originating at (y = 0, t' = t) and ending at (y = 0, 
t* = 0), after reversing the direction of time.) 

Introducing a turbulent diffusivity into equation (III-16) yields the 
trajectory 

Y(f) y+ v^(w^, 


For this trajectory 


X(t) 


0 


t 


Y(t')dt' 


has mean 


and variance 



nyt 

= 


(X - X)^ = 



(Wf -f W^) -f w,) df dt 

min(t',t") - df dt" = 

t 0 
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This mean and variance agree with those of equation (III-15), showing process 
(III-16) to be the correct form of a Wiener process with two fixed end 
points. 

There is an advantage to approaching the calculation of conditional 
moments via prescription of this modified Wiener process. The stochastic 
process (III-16) contains more information than the pdf (eq. (III-15)), for 
the latter describes the stochastic process only at a single time. The 
joint distributions at multiple times can be inferred from the process (eq. 
(III-16)). Or more generally, its description as a continuous dynamical 

process makes W^, amenable to analysis. 

Reflecting boundary . - The case where a reflecting boundary is inserted 
at y = 0 can also be dealt with. This has already been discussed in chap- 
ter II (eq. (11-37)). When boundaries are present in turbulent flow, they 
make the turbulence inhomogeneous; so it is not entirely consistent to take 
K constant and apply a reflection boundary condition. However, for illus- 
trative purposes it is an interesting case to consider. 

As mentioned in chapter II the reflection condition is incorporated by 
taking an absolute value. Thus the stochastic process 


Y(t) 




(III-18) 


satisfies a reflection condition at y = 0. For convenience in equation 
(III-18) the source yg of the trajectories will be put at the boundary 
yo = 0. Then equation (11-37) gives the pdf for the process (III-18) as 


P(y,t) 


From this the moments 


1 g-y^/4Kt 
TrKt 


H(y) 


y = 2 y'kt/iv 



(III-19) 


(III-20) 


are computed as in equation (1-4). 

By equations (III-18) and (III-6) 

X(t) = n y^^^(W^,|df (III-21) 

Because IW^l is not Gaussian distributed, X is also no longer Gaussian; but 
its moments can be computed from the known properties of W^. For example. 
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J = n 



dt' 


= 2n 



df 




(III-22) 


2 

and X is given by 


= 2n^K IW. ,W.„| dt" df = 4n^K f f^ |W^..W.„1 dt" df (III-: 

I t t *'0 *'f t t I 

To evaluate the average occurring inside this double integral, the pdf of 
W^i and Wfi for t" 2 is required. Since W^i and Wf are 


t y. t 


23) 


jointly Gaussian with , = f , = t", and W^,W^„ = f , this pdf is 


P ( I » II ) — 


^^p|'2f (t" - f ) [^f t" (^t" 


2tt yf (t" - f ) 

(eq. (5.1.28) of ref. 8). From equation (III-24) and the definition 


(III-24) 


■^f *^t " I “ SS I '^t ' '^t " I ^ ^ \ ’ ”t " ^ '^'^f 


dw. 


f‘ 


it is found that 


^ 2f ^^^(t" - f )^^^ 


irt" 


. sin (cos ^ Vt'/t") - 

Vf /t" cos ^ 

3 / ^ 

sin (cos 

VtVtM 


which is a comparatively complex formula. Substituting this into equation 
(III-23) and changing variables of integration to e and <p where 


give 


cos 6 = ^t ' / 1 " and cos <p = -^f /t 
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= 3n^Kt^/4 

The corresponding variance is 


(III-26) 



Clearly, the reflection condition could also be imposed on the process 
with two fixed end points (eq. (III-16)). In this case 


Y(f) 




or, evaluating this on y = 0 for simplicity. 


Y(f) = 



(III-27) 


Corresponding to the stochastic process (III-27), 


X(t) = n 



t 




(III-28) 


Since W^, - W^t'/t is a Gaussian random variable with variance 

t'(l - t'/t) (eq. (III-17)), the mean value of jw^, - W^t'/t| is computed to 

be ^2t ' ( 1 - t ' /t) /tt. Hence 


X = 2n 



(III-29) 


"T -2 

This result, as well as the value of the variance X - X , is given in ref- 
erence 1 (eq. (5.103)). 


Example 2: The Surface Layer 

Criterion for eddy diffusion . - Generally, the condition for turbulent 
dispersion to be described by a diffusion process is that the correlation 
time T|_ of turbulent velocity fluctuations can be treated as being negli- 
gible. In most cases this is possible at times that are long compared with 
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T|_, but in the early stages of dispersion this will only be possible if Tl 
is negligible at the position of the contaminant source. 

Now, a vanishing correlation time is associated with a vanishing corre- 
lation length of turbulent eddies. In turn, is determined by a 
length scale of the geometry within which the turbulence exists. Near a 
plane wall, since distant regions of the flow do not affect the near-wall 
turbulence, the only geometric scale is distance from the wall. Thus the 
correlation length here is proportional to distance from the wall; in parti- 
cular it tends linearly to zero at the wall. This near-surface region, in 
which “ y, with the wall at y = 0, is called the surface layer . 

The conclusion to be drawn from the two paragraphs above is that eddy 
diffusion is valid at short times only for sources located at the wall: the 

argument being that the diffusion approximation is valid when t >> Ti_, 
and because T[_ 0 at the wall, this criterion is satisfied by surface 

sources for all time. The effects of finite correlation time (or, the case 
of elevated sources) are discussed in chapter IV. 

The surface layer . - The surface layer will be analyzed as a semi- 
infinite field of turbulence above a plane wall, in which the eddy diffusiv- 
ity increases linearly with height. This form of diffusivity follows from 
dimensional analysis, for the only velocity scale for wall turbulence is the 
surface friction velocity u* and the only length scale is y. Thus 

K(y) = eu*y (III-30) 

where the constant of proportionality e is usually written as 
e = 0.4/(Pr)j, (Pr)j being the turbulent Prandtl number (ref. 2, p. 51). 

Our analysis is simplified by introducing as the independent variable 

t = Bu*t, where t is actual time and t has the dimension of length. Then 
equation (III-l) with diffusivity (III-30) becomes 


dY = dt + -^/2Y dW^ (III-31) 

and, for a surface source, Y(0) = 0. By equation (III-31) and Ito's theorem 
dY^ = nY^"^ dt + n(n - 1 )y""^ dt + nY'^"^ dW^ 
which, after averaging, gives 

dY'^ = n^ y"^”^ dt (III-32) 

"h h 

for the n*"" moment of Y. By induction on n, the solution to equation 
(III-32) satisfying y'^(O) =0 is 


y" = r(n + l)t'^ (III-33) 

where r is the gamma function: r(n + 1) = nJ when n is an integer. 


37 



A knowledge of all of the moments of a random variable is as good as a 
knowledge of its pdf; in principle one determines the other. Because 


r(n 



x"e-" 


dx 


it is clear by inspection that the pdf having the moments (III-33) is 


P(y,t) = |e H(y) 


(1 11-34) 


(As previously, H(y) is the Heaviside unit function.) This is the one and 
only probability density function, all of whose moments are given by equation 
(111-33). 

The general method for determining a pdf from its moments is to note 
that the Fourier transform of P(y,t) is (if the sum converges) 

00 

P(y,t)dy v"(t) (III-35) 

n=0 



so, by inverting the transform, 

=57/" At)dk 


n=0 


(III-36) 


In the case of equation (III-33) 



(ik)*^ - 1 

n I (1 - ikt) 


which has the pdf (III-34) as its inverse transform. Formula (III-36) is 
used later in our analysis. 

The streamwise mean velocity of a turbulent boundary layer near a plane 
wall has a logarithmic form 


U(y) = 0.4 u* In 



(III-37) 


where the constant yg is the effective surface roughness height. Equation 
(III-31) follows, again, from dimensional arguments (ref. 2, p. 155). As in 

the last example, we will next calculate X and X from equation (III-2). 
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Since dX = U(Y) dt. 


X = UTyT dt = (Pr)^ Y - In dt (III-38) 

The value of In Y is easily determined. Equation (III-33) is valid for all 
real values of n, so by differentiating it with respect to n and setting 
n = 0 


dn 


n = 


0 


In Y = In t 


Y 


(III-39) 


where 


Y = 0.5772156... 


dr(x) 

dx 


x=l 


is Euler's constant. Thus (for convenience, from here on we set (Pr)j = 1). 


X = 



(III-40) 


The variance of X is 


- X^ 




In Y(t') In Y(t") - In Y(t') In Y{t*'T dt" dt' 


t 


The integrand of equation (III-41) can be evaluated analogously to equation 
(III-39); an expression for 


f In Y(t') In Y(t") - In Y(t') In Y(t") dt" dt' (III-41) 



fnm(t',t") = Y"(t') Y"^(t") - Y"(t') Y"’(t") (III-42) 


is derived and differentiated with respect to n and m, setting n = m = 0 
in the result. Details of this procedure are given in appendix A of this 
chapter. The result 
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(III-43) 


was derived by reference 11 through an analysis of the Fokker-Planck (i.e., 
diffusion) equation (III-4). 
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Conditional mean of X . - As in the previous example the mean value of 
X at fixed y can be determined (cf. eq. (III-15)). Recall that this mean 
is calculated from P(xly): 


X(y,t) 



P(xly,t)x dx 


(III-44) 


Because X(t) = 1 n[Y(t ' )/yQ]dt ' , equation (III-44) can alternati vely be 


written 


X(y,t) = 


//■ 

*^0 


P(y‘ ,t'| y,t)ln dy' dt ' 


It is convenient to use the fact that P(y‘ly)P(y) = P(y',y) and calculate 
[X(y,t) - X(t) ]P(y,t) , where X(t) is as in equation (III-38). This quantity 
is given by 


[X(y,t) - X(t)]P(y,t) 


-ff 


P(y‘ .t' ,y,t)(ln y' - In t + Y)dy‘ dt' 


(III-45) 


But, by expanding the y-dependence as in equation (III-36), the right side 
of this last equation becomes 


\_ 
2 TT 



g-iky\ ^ ( . ik)* ! Y^(t) In Y(t') - Y^(t) In Y(t')dk dt ' 
n=0 


/•t /•“ _ “ ■ 

hj 

•'O n=0 


^ ^ ^). — f ( t t ' ) 
nl dm ^nm^**’^ ^ 


dk dt' 


m=0 


(iri-46) 


with f|^|^ defined in equation (III-42) and given explicitly in appendix A 
by equation (III-A7). In appendix B equation (III-46) is evaluated to show 
that 


X(y,t) - X(t) = t[e^E^(0 + In C + y - 1] (III-47) 

CD 

where K = y/t and E, (5 ) = / e~^/x dx is the exponential integral. Ex- 

pression (III-47) was derived by Chatwin (ref. 11). His approach is simpler 
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than that used here, but our purpose is to illustrate the analysis of sto- 
chastic processes without resort to the diffusion equation (which is not 
available for non-Markov processes). For Markov processes the diffusion 
equation is available, so analysis can be based either on it or on the cor- 
responding stochastic differential equation. The present example shows the 
equivalence of these analyses. In a sense this can be called an equivalence 
of the Lagrangian and Eulerian approaches to turbulent transport. 

Streamwise dispersion . - In most experiments the surface layer consists 
of the near-wall region of a turbulent boundary layer. In these experiments 
the source of contaminant is steady and dispersion takes place in the cross- 
stream direction y as a function of streamwise distance x from the source 
at X = 0. Thus these experiments should be modeled by a stochastic process 
Y(x) with X the independent variable. When dt = dx/U(y) is used to change 
from t- to x-dependence, equation (III-31) becomes 


= U[Y'(x)] |^U[Y(x!] "^^x 


(III-48) 


The corresponding Fokker-Pl anck equation is (cf. eq. (11-16) with x re- 
placing t) 


3P(x,y) ^ i_ y i_ P(x,y) 
3x ay ^ ay U(y) 


(III-49) 


One usually defines a normalized concentration C(x,y) by P(x,y) = 
U(y)C(x,y), which shows the probabilistic meaning of C(x,y). The concen- 
tration C(x,y) satisfies 


Ufvl = i- V 

ax ay ^ ay 


(III-50) 


If U(y) is given by a power law, U(y) “ yP, then a change of indepen- 
dent variable to YP‘'’1 (x) transforms equation (III-48) to the form of equa- 
tion (III-31), and it can be solved as above. However, when the more correct 
logarithmic form (III-37) is used for U(y), equations (III-48) and (III-49) 
can no longer be solved in closed form. In this case a finite-difference 
version of equation (III-48) can be solved with relatively little computa- 
tional effort by the methods described in chapter II. Reference 12 numeri- 
cally solves an equation like equation (III-48) downstream of a line source 
located on the wall below a turbulent boundary layer. From an ensemble of 
solutions, statistical moments of Y(x) are evaluated. Data on plume moments 
downstream of a source in a boundary layer are provided in reference 13. The 
plume moment is defined as 


Yp(x) 


= ^ yC(x,y)dy = y ^P(x.y)dy 


(III-51) 


Numerical results of reference 12 for Yp(x) are shown in figure III-2, along 
with measurements of Yp(x). It can be seen that the eddy diffusion model 
does a good job of describing the evolution of Yp. 
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Figure III-2. - Increase of plume moment with dis- 
tance downstream of wall source in turbulent 
boundary layer. Distances are normalized by 
boundary layer thickness 6. 


Example 3: Channel Flow 

We analyse this example in considerably less detail than the previous 

two. 

Consider the turbulent flow through a channel with walls at y = 0 and 
at y = 1. A diffusivity suitable to this flow is 

K(y) = u*6y(l - y) (III-52) 

with a as in example 2. Again if we let t = u*et, Y(t) satisfies 


dY = (1 - 2Y) dt + y2Y(l - Y)dW^ (III-53) 

The Fokker-Planck equation corresponding to equation (III-53) is 


af _ 

at " 


8y 


[y(i - y) |y_ 


(III-54) 


For the reasons given at the beginning of the previous example the eddy dif- 
fusion model, equation (III-53) or (III-54), of dispersion is only suitable 
for wall sources. Hence, we consider the initial condition Y(0) = 1 to 
equation (III-53). The corresponding condition to equation (III-54) is 
P(y»0) = <s(l - y) » where 6(*) is the Dirac delta function. 

By Ito's theorem 


dY" 
dt ■ 


-n(n + 1 )y" + 




(III-55) 


subject to y" (0) = 1, determines the moments of Y. These moments are 
readi ly seen to be 
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Y = I (1 + e"^*) 


-2t -6t 


= i + ^T5 — + - 


(III-56) 


/„i \2 X 2m + 1 _-m(m+l)t 

^ 2^ n- - mi n' + m + IJ ® 

m=0 


That Y*^ satisfies equation (III-55) is seen by substitution; that the 
initial condition is satisfied is seen by writing Y^ (0) as 


(n.*r 


(n + m + 1 


n + m + 



m=0 



-E 

m=l 


n + ml n - ml 


= 1 


The expression for Y(t) is shown in figure III-3 with data from Sullivan 
(ref. 14). He measured the trajectories of neutrally bouyant particles in 

turbulent channel flow, from which Y(t) was determined directly. 



Figure III-3, - Mean trajectory in turbulent channel 
flow of particles released at wall. 
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The probability density for this example is obtained most simply by 
solving equation (III-54): 


1 1 

P(y,t) = ^ - l)(2n + 1) (III-57) 


n=0 


where ^n(* ) the n^h Legendre polynomial. The moments (III-56) 
follow from the pdf (III-57), showing the equivalence of the moment and pdf 
representations of the statistical properties of Y. 
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APPENDIX A: EVALUATION OF X-VARIANCE FOR EXAMPLE 2 


The procedure for determining the variance of X consists first of 
finding an expression for fnm with n and m real, then differentiating 
with respect to n and m (cf. eq. (III-41) et seq.)* Because of the non- 
anticipating property of Y(t), Y(t") is uncorrelated with dWt< when 
t" < t'. Therefore multiplying 

dY'^(t') = n^Y'^“^(t')dt' + nY"“^(f )dW^, (III-Al) 


by Y'’i(t") and averaging give 


^ Y"(f)Y^(t") = n^ Y""^(t')Y"(t") 


(III-A2) 


where t" is treated as constant and t' > t". Since Y^(t') satisfies 
equation (III-32) with respect to t', fnfj,(t',t") as defined in equation 
(III-42) satisfies 


df ^(t',f) , 

QID n^f ft' t"l 


(III-A3) 


Because t' > t" in equation (III-A3), f^m must satisfy an initial condi- 
tion at t' = t"; clearly 


f^^(t",t") = Y'^‘^'^(t") - Y'^(t") Y"’(t") (III-A4) 


with the moments of Y given by equation (III-33), is this condition. Also, 
the recurrence relation (III-A3) starts from 


= 0 


(III-A5) 


When n and m are integers, it can be seen that 

n 


nm ^ 


p=0 


p.'^n - p.' 


(m + pi - mip.' ) 


(III-A6) 


satisfies equations (III-A3), (III-A4), and (III-A5). 

Expression (III-A6) can be extended to noninteger n and m by re- 
expressing it as a contour integral: 


f 


nm 


r^( n + 1 ) 
2iri 



(f)m-S(f _ t")^~"^(-l)^r(s) 
r(n + 1 + s)r(l - s) 

x [r(m - s + 1) - r(m + l)r(l - s)]ds 


(III-A7) 
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The contour of integration envelopes the poles on the negative real axis as 
shown below. 


Original path 


Complex x-plane 




Distorted path 


(This contour is appropriate when t" < t'/2. We appeal to analytic continu- 
ation to make equation (III-AIO) valid for all t" < t'.) 

As described in the text below equation (III-71), 


-2 

= 2 


Differentiating equation (III-A7) gives 


d^f 


nm 


dn dm 


r 

r 

d^f 

nm 

Ij 


dn dm 

0 gives 


r(s)r 

■(1 + 


dt" df 


n=m=0 


(III-A8) 


-hf CMi- 


s) + v]ds (III-A9) 


where r'(x) = dr(x)/dx, vp(x) = r'(x)/r(x) and -y- = -'l'(l). Upon evaluating 
the residues of equation (III-A9), the right side becomes 


OO OO p 

/ f \ P Ml + p) + Y V 1 

p " prVt"-t7 

p=l p=l r=l 


r=l p=r 



t7(f-t") 


ln(l + x) 
x(l + x) ■ 


dx 


(III-AIO) 


where use has been made of the identity 
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+(1 + p) 


p 



for integer p. Equation (III-AIO) is the integrand of equation (III-A8), 
which therefore becomes 


/*^ /* ^"/(t'-t") , ,, ^ \ 

1 1 1 


0 *^0 


= 2t‘ 




where t = t'/t, a = (f - t")/t. This triple integral is quite harmless: 
after integrating by parts with respect to a it becomes 



ln(T/A) 
T / a - 1 


da dr 



X In X 

1 - X 


dx 


/ I X 2 

In X dx - / dx = ^ - 1 (III-A12) 

u *'o ^ ^ ® 


The substitution x = a/t was made in the first step. Equation 
(III-43) is thus obtained. 
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APPENDIX B: EVALUATION OF 7’(y,t) FOR EXAMPLE 2 


We wish to evaluate expression (III-46) with given by equation 

(III-A7). Differentiating equation (III-A7) with respect to m and setting 
m = 0 give 


dm 


r^(n 


m=0 


2iri 


11 


I 


u- 




n+s 








[♦(1 - s) + Y]ds 


■ ♦ P) ^ r] (III-Bl) 

P=1 


By substituting this and exchanging orders of integration, equation (III-46) 
becomes 




n 

nJ(ik)" ^ 


p=l 


’t'd + P) + Y 
pi n - pi 



(t - t')'^ P t'P dt' 


dk 


(III-B2) 

Because the last integral is equal to (n - p)lpl/(n + 1)1, this ex- 

pression equals 


:L f 

2ir J 


but 


n=l p=l 


n p n n 


(III-B3) 


4»(p + 1) + Y = i . (n + l)(*(n) + t) - n 

p=l p=l r=l r=l p=r 


so expression (III-B3) equals 

cx> 

n=l 




dk 



e-xy/t 


ln(l - x) 
x(l - x) 



dx 


(III-B4) 
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where x = ikt. The residue for the 1/(1 - x) in square brackets is just 

e~^ (III-B5) 

where 5 = y/t. The first term can be evaluated by distorting the path of 
integration to envelope the part of the positive real axis greater than 1, 
as shown in the following sketch: 


Imag (s) 





Complex s-plane 


Real (s) 


Then this term becomes 



^x ln(l - x) 
^ x(l - X) 


dx 


Since In(-x) = ln|x| - ni above the real axis and In(-x) = ln|x| + Tri be- 
low the real axis, this last integral equals 



(III-B6) 


provided the integrals are interpreted in the generalized sense (ref. 15); 
that is. 




In X dx = - In 5 


- Y 


(III-B7) 


The first integral in equation (III-B6) is just the exponential integral 
Ei( 5). Thus, collecting the results (III-B5), (III-B6), and (III-B7) gives 

e'^ (In c + Y - 1) + El( ?) (III-B8) 


for the value of the expression (III-46). As given in equation (III-34), 
P(y,t) is e“^/t and expression (III-B8) is the right side of equation 
(II 1-45). Thus we obiain equation (III-47). 
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INTERLUDE: A POINT CONCERNING MODELING 


In examples 2 and 3 of chapter III it is noted that the eddy diffusion 
model can only be used for dispersion from surface sources. For elevated 
sources this model is incorrect in principle and in practice. It is some- 
times suggested that a time-dependent diffusivity could be used to rectify 
this failure (ref. 1, section 5.8). In this procedure, time must be referred 
to an origin at the contaminant source. However, the justification for this 
suggestion is usually based on a false identification of diffusion processes 
and stochastic processes with Gaussian pdf's. The Gaussian distribution 
solves a diffusion equation with time-dependent diffusivity, as a matter of 
fact, but this does not make Gaussian processes diffusion processes. In non- 
homogeneous turbulence the pdf is not Gaussian, so the failure of eddy dif- 
fusion there certainly cannot be rectified by a time-dependent diffusivity. 

An illuminating discussion of the irrationality of time-dependent dif- 
fusivities is given by Taylor (ref. 16). He remarks that a diffusivity 
which is a function of time violates the principle of superposabi 1 ity: if 

two sources of contaminant start at different times, then where their clouds 
overlap it becomes necessary to assume that the diffusivity has two values, 
each referred to a different origin at the appropriate source time. This is 
nonsensical since the material emitted by each source is indistinguishable 
from that emitted by the other. 

However, this scheme for dealing with elevated sources is illogical for 
another reason, which has to do with a general principle of turbulence model- 
ing. Namely, the coefficients in a dispersion model can only be functions 
of the turbulence and not of the contaminant distribution (which, after all, 
is what is to be predicted.) Hence in stationary turbulence these coeffi- 
cients cannot be functions of time; for example, in chapter III it would be 
incorrect to let K(y) be a function of time. Even when the turbulence is 
not stationary, the coefficients should be referred to the time origin of 
the turbulence and not to any source time. The next chapter describes a 
model of dispersion from elevated sources in stationary turbulence that 
abides by this principle of modeling. 
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CHAPTER IV: DISPERSION WITH FINITE TIME SCALE 


In chapter III examples in which particle trajectories are modeled by 
Markov processes were considered. In such processes the random component of 
the velocity is represented by white noise (recall that the white-noise 
process is dWt/dt). Because white noise has a zero correlation time 
scale (cf. eq. (II-5)), by representing the randomness in this way one is 
assuming that the correlation time of velocity fluctuations T|_ can be 
regarded as negligible. In examples 2 and 3 of chapter III it was noted 
that the Markovian eddy diffusion model was appropriate for dispersion from 
surface sources. The reason is that in the nonhomogeneous turbulence above 
a wall the velocity time scale tends to zero as the surface is approached. 
However, away from the surface this time scale is not negligible: dispersion 

from elevated sources must be modeled by a process with a finite correlation 
time for velocity fluctuations. 

At high Reynolds numbers (Re = u'®'/v) the time scale of acceleration 
fluctuations is of order (Re)~^'^. This is so because acceleration has 
its rms value determined primarily by those components of turbulence that 
have a short time scale compared with the velocity time scale T[_. Analysis 
in reference 5, section 18.3, shows that the acceleration time scale T^ 
i s of order (Re)~^/^TL. Thus, while Ti_ is of order 1, T^ tends 
to 0 like (Re)“l/2 at high Reynolds numbers. Consequently it is suitable 
to model acceleration by a white-noise process; indeed, in writing equation 
(I-l), terms of order (Re)“^'^ were ignored, so letting Tg ~ 0 is 
consistent with that equation. 

An equation in which acceleration is modeled by white noise is the 
Langevin equation (ref. 7): 


dv(t) = dW^ (IV_i) 

Here _ and T|_ are constants and V(t) is velocity. The initial 
condition for equation (IV-1) is that V(0) is a Gaussian random variable 

2 

with mean zero and variance a^. The Langevin equation has its origin in the 

theory of Brownian motion, but it was used in a discretized form by Taylor 
(ref. 17) as a model of turbulent dispersion. Application of the analysis 
in this chapter to environmental dispersion can be made along the lines of 
reference 1, chapters III and V. 


Langevin Model for Homogeneous Turbulence 

Equation (IV-1) is a Markovian differential equation of the type 
discussed in chapter II, except that now it determines particle velocities 

V(t). The trajectory Y(t) = / V(t') dt' is not a Markov process 

*'0 
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(although V and Y jointly are a Markov process in phase space, ref. 3). 
Because the turbulent velocity has a finite time scale, particle trajectories 
must be nonMarkovian. But, by retaining the Markovian property at a higher 
level - that is, by making the velocity be a Markov process - the dispersion 
model remains tractable. There is no general theory for completely nonMarkov 
stochastic processes, so it is useful (for numerical as well as theoretical 
analysis) to maintain this connection to Markov processes at some level of 
modeling. 

Equation (IV-1) is a linear equation for V(t). Because dW^ is 
Gaussian and equation (IV-1) is linear, V(t) is also Gaussian; thus V(t) is 
characterized by its mean and its variance. Averaging equation (IV-1) gives 

dV -7 


But V(0) = 0, so the mean of V(t) is zero. By Ito's theorem 


1 ^ 
2 dt 



T 


L 


Because V(0) is Gaussian with variance 


a 


2 

V’ 



(IV-2) 


is the solution of this equation for the variance V(t). Thus V(t) is 
Gaussian with constant mean and variance; in fact, it is a statistically 
stationary random process. The pdf of v is 


P(v) 



(IV-3) 


In stationary, homogeneous turbulence the particle velocity is statisti- 
cally independent of particle position; hence it is a stationary random 
process. For this reason equation (IV-1) is a model of dispersion in homo- 
geneous turbulence: the constancy of and T[_ makes it a model of 
homogeneous turbulence. 

Carrying the analysis of second moments of V one step further, we 

2 

calculate the correlation function R(t) = V(t + t )V(t)/a^, t ^ 0. With t 
fixed and t variable, V(t + t) satisfies (cf. eq. (IV-1)) 

dv(t + t) = dT + ^ dW^^ (IV- 4 ) 
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By the nonanticipating property, V(t)dW. . = 0, t ^ 0. Hence multiplying 

equation (IV-4) by V(t), dividing by a^, and averaging give 

dR(-c) _ -R(t) 

dT T! 


This has the solution 


R(t) = e 


-WT, 


(IV-5) 


which satisfies the initial condition R(0) = 1. Thus the role of Ti_ as 
a correlation time scale is clear. Here it could be an integral time scale. 



R(t )dr 


but it is more natural to regard it as a scale for instantaneous decorrela- 
tion rate: 


T 


L 


1 

WT 


dRllI 

dt 


-1 


1=0 


(IV-6) 


The reason for the interpretation (IV-6) of T[_ is that when the turbu- 
lence is stationary but nonhomogeneous, V(t) is not a stationary process; in 
this case an integral time scale is neither clearly defined nor a property 
of the turbulence. (Remember that V(t) is the velocity of a particle. 

When the turbulence is not homogeneous, this will be statistically dependent 
on particle position and consequently will not be stationary. A time- 
dependent integral scale could be defined, but it would not be a property of 
the stationary turbulence.) Equation (IV-6) defines an instantaneous prop- 
erty of V ( t ) ; therefore in nonhomogeneous turbulence it defines a local 
property of the fluid flow. If we wish to use equation (IV-1) for nonhomo- 
geneous turbulence, equation (IV-6) must be used to define a local T[_(y); 
similarly ay(y) would then be related to the local value of the turbulent 
velocity variance. Later the application of the Langevin equation to non- 
homogeneous turbulence will be considered; first we explore the case of homo- 
geneous turbulence further. As in example I of chapter III, we consider dis- 
persion in a uniform shear flow but now allow for the finite time scale of 
the turbulence. 


Shear Dispersion 

The random trajectories of fluid particles are given by 

dY(t) = V(t)dt 
dX(t) = nY(t)dt 



(IV-7) 
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(cf. eq. (III-6)) with V(t) determined by equation (IV-1) and the mean 
x-velocity being U(y) = ny. The turbulent velocity in the x-direction is 
ignored with respect to the mean; it could easily be included as a random 
process analogous to, but independent of V(t) - or, to model Reynolds 
stresses, it could be correlated with V(t). 

Because V(t) is Gaussian and equations (IV-7) are linear, X(t) and 

Y(t) are Gaussian. Furthermore, if they are subject to the initial condi- 
tions X(0) = Y(0) = 0, then "X = T = 0. We proceed to analyze the second 

moments of X and Y. 

Integrating the first equation (IV-7) with Y(0) = 0 gives 



I 


SO 



V(f)V(t")dt" dt' 


= 2 



V(f )V(t")dt" dt 


(IV-8) 


By equation (IV-5) V(t' )V(t") = exp[-(t‘ - t")/T|^], for t' > t". After 
this is substituted, the integrals in equation (IV-8) can be evaluated: 


y2(t) 


, 2,2 

2»v\ 




(IV-9) 


As has been said previously, when t >> T[_, the Markovian models of 
chapter III ought to be valid. When t >> T[_, equation (IV-9) reduces to 





t 


2 

Comparing this to the Markovian example 1 of chapter III, where Y = 2Kt, we 
find that 


K = a^T|_ (IV-10) 

is the eddy diffusivity of the Markov process to which equations (IV-1) and 
(IV-7) reduce as t/lL “. Because Y(t) is Gaussian, the fact that its 
variance tends to that of example 1 of chapter III shows that it 

asymptotically has the Markovian pdf. At short times, Y^ + t^. 

The variance of X is 
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dt" 


(IV-11) 



Y(t' )Y(t")dt' 


in which 


^11 

Y(t')Y(t") = / f V(x')V(x")dx' dx" 

*'^0 •'o 


= a 


2 f ^" -|x'-x"|/Tl 


dx' dx" 


- »vT|_ I e + e - e - 1 * T 


(III-12) 


and t" £t'. Substituting equation (IV-12) into equation (IV-11) and 
integrating give 



9 2 2,4 
2n 





(IV-13) 


As t/T -► “ this reduces to equation (III-IO) with K given by equation 

(IV-10). As t/T^ + 0, a^n^t^/4, so initially 7^ grows quite 

~~2 ~~2 
slowly compared with Y , although ultimately it grows more rapidly than Y . 

The covariance XY is readily computed: 

XT = n Y(t)Y(t')df 

*^0 



- 1 + e 


-t/T, 


(IV-14) 


having used equation (IV-12). Because X and Y are jointly Gaussian 


P(x,y) 


exp 

2 2 2 2 — 
-x'^ Y'^ _ y^ + 2xyXY 

2(X^ Y^ - XT^) 


2tt - XT^ 


(IV-15) 
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(ref. 8, eq. (5.1.28)) with X , Y , and XY given by equations (IV-9), 
(IV-13), and (IV-14). The function P can be regarded as the concentration 
of a contaminant released initially by a point source at x = y = 0. The 
numerator of the exponential in equation (IV-15) is constant on ellipses in 
the x,y-plane: these ellipses are the surfaces of constant concentration. 

At short times their major axis is nearly along the y-axis; at long times it 
is nearly along the x-axis. 

The pdf (IV-15) corresponds to the pdf (III-13), and 


P(x|y) = P(x,y)/P(y) 



(IV-16) 


corresponds to equation (III-15). Equation (IV-16) is the probability density 
for the x-position of those trajectories that at time t are at a given 
y-position. It is also the normalized streamwise distribution, at fixed 
cross-stream position, of contaminant released by a point source. 

From equation (IV-16) the first and second moments of X(t) for those 
trajectories that have the fixed end point Y(t) = y are 


X(y,t) = y ^ = 

2 

v/t. 


exactly as was true for equation (III-15), and 


(IV-17) 


X^(y,t) - X^(y,t) 



? 2 2t-4 
= 2n 





(IV-18) 


2 —2 

As previously, X (y,t) - X (y,t) is independent of the end point y. When 

2 2 5 

t/T| + 0, expression (IV-18) tends to n a t /60T, , which is much smaller 

L 2 2 4^^ 

than the unconditional variance n a^t /4 (cf. eq. (IV-13)). The reason for 

the conditional variance being smaller than the unconditional variance is 
that the conditional variance is the spread of contaminant about the mean 
(eq. (IV-17)), while the unconditional variance is the spread of these 
spreads. 
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The preceding analysis of dispersion by a linear shear in homogeneous 
turbulence describes horizontal dispersion in a boundary layer. As mentioned 
in chapter III, this model is not suitable for vertical dispersion near a 
boundary because the turbulence is not homogeneous in that direction. Thus, 
while a reflection boundary condition could be imposed, following the ap- 
proach of chapter III, example 1, such analysis is not strictly justified 
and will not be pursued here. Instead we proceed directly to a discussion 
of dispersion with finite time scale in a nonhomogeneous medium. 


Dispersion in Nonhomogeneous Turbulence 

In nonhomogeneous turbulence the coefficients in equation (IV-1) must 
be made functions of Y. For the present consider the case with T|_ a 
function of y and ay constant; later we will consider the case in 
which ay is also a function of y. 

Now equation (IV-1) has the form 


dV 


V 

"T^ 


dt + 



(IV-19) 


Along with dY = V dt the model (IV-19) describes trajectories through non- 
homogeneous turbulence. Equation (IV-19) is to be solved subject to the 
conditions that Y(0) = h and that V(0) is Gaussian with mean zero and 
variance ay for a source located at y = h. (The initial condition on 
V follows from the requirement that contaminant particles move with the 
fluid. At the point of release the fluid velocity is assumed by equation 
(IV-19) to be Gaussian with mean zero and variance ay.) Generally equa- 
tion (IV-19) must be solved numerically (refs. 12 and 18); however, here it 
will be solved by asymptotic analysis. 

It is obvious that at sufficiently short times T[_ can be considered 
to be constant with value T[_(h). During this initial period, dispersion 
takes place as described previously for homogeneous turbulence. To delimit 
and go beyond the initial stage, one might construct a solution to equation 
(IV-19) as a power series in t and W-t. However, a different, though 
similar, approach will be followed here: we consider dispersion in a slowly 
varying medium (ref. 19); that is, we assume Tl can be written as a 
function of ey where e « 1. Our objective is to determine how inhomo- 
geneity affects the solution to equation (IV-19). 

With the functional form Ti_ (ey), equation (IV-19) reduces to the 
equation for a homogeneous medium when e = 0. For small but nonzero e 
the solution to equation (IV-19) can be expanded in a power series 


V(t) = Vg(t) + eV^(t) 
Y(t) = YQ(t) + eY^(t) 


(IV-20) 


The terms with subscript 1 (and greater) are due to the inhomogeneity of the 
turbulence. 

The equations solved by Vq and Vx are found by substituting the 
expansion (IV-20) along with the Taylor series 


57 



(IV-21) 


TjeY) = 1^(0) + eTL(0)YQ . . . 


where T|^ = dT|^/dy, into equation (IV-19). Without loss of generality the 

source height h has been set to zero. Equation (IV-19) then can be sepa- 
rated into a set of equations by equating to zero the coefficient of each 
power of G. To order g this results in 


and 


“''o = - T[f<5T 


•Vi 


T~ 
W)° 




dVl = 


''l .. + ^L(O) V u Ht 

\loy Vo 




(IV-22) 


(IV-23) 


where dYg = Vq dt and dY^ = Vj dt. Clearly the solution to equation 
(IV-22) is that given in the previous section; in particular 




9 2,2 
Wi 





> 


J 


(IV-24) 


where it is understood that Ti_ is evaluated at y = 0. 

At order g the mean value of Y is no longer zero. Upon averaging, 
equation (IV-23) becomes 



^1 + 1 ^^0 
T, 2 -2 dt 

^ 'l 


1 2 

In this equation YqVq has been written in the equivalent form dYg/dt. 

Substituting equation (IV-24) and solving this, subject to V^(0) = 0, gives 



2t' 



Correspondingly 







-t/T^ 



(IV-25) 
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Thus the effect of inhomogeneity is to cause a nonzero mean vertical veloc- 
ity. In physical terms the source of this mean velocity is that particles 
moving in a direction of increasing T|_ maintain their direction of motion 
longer on average than particles moving in a direction of decreasing Tl; 
thus particles drift up gradients of T. . 

L 2 ' 

As t/Tj^ the velocity (IV-25) becomes If K is given by 

equation (IV-10) with Tl a function of y, this velocity is 


dt dy 


(IV-26) 


In the eddy diffusion model, equations (III-l) and (11-19), dY/dt is also 

dK/dy; therefore as tUi -*■ equation (IV-25) again reduces to the eddy 
diffusion result. In the next section it will be shown more completely that 

the stochastic process Y = Yq + eY^ tends asymptotically to an eddy dif- 
fusion process. 

Although dY/dt tends asymptotically to the formula (IV-26), it must 

initially be zero, for Y(0) = 0. The effect of the finite correlation time 

scale is that particles "remember" that they had this zero initial velocity, 
and their mean velocity only gradually becomes nonzero, tending eventually 
to the velocity (IV-26). In the eddy diffusion model the memory time scale 
is zero, so random trajectories immediately obtain the mean velocity (IV-26). 



t/iL 


Figure IV-1. - Comparison of mean rise of particles as 
given by equation (IV-25) with that given by eddy dif- 
fusion model. 



Figure IV-2. - Behavior of plume mo- 
ment with distance downstream of 
elevated source in turbulent bound- 
ary layer. Distances are normalized 
by boundary layer thickness 6. 
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This is illustrated in figure IV-1 by a comparison of the expression for Y, 

2 I J- 

with the diffusion result Y = OyT|_t. 

Shown in figure IV-2 is a numerical solution given in reference 12 of a 

model like equation (IV-19) for Y(x) (actually, for Yp(x), see eq. (III-51) 
et seq.) in a turbulent boundary layer; for comparison,^data from reference 
14 are included. Here downstream distance x plays the role of t in the 

preceding. In both theory and experiment, dY/dx initially is zero and Y(x) 
only gradually becomes parallel to the eddy diffusion asymptote shown by the 
dashed line. 

“7 

Next consider the mean squared value of Y: it will be shown that Y 
is unaltered to order e. Since Y^ = Yq + eYj^Yq + O(e^) and 



we need only show that V^Tt^TV^Tt^’T = 0. When t' > t", multiply equation 
( IV-22) by Vj(t") to find 

dVQ(f)V^(f) Vo{'f)V^(t") 

_ = ^ 


and when t‘ < t", multiply equation (IV-23) by VQ(t') to find 


dVort')Vi(t") 

dP 


VQ(t')Vj(t") T|^ dY2(t")V„(t') 

’I 


? 2 

That YQ(t")VQ(t') s 0 is clear from symmetry: Yq is symmetric with 

respect to reversal of the sign of in equation (IV-22), while Vq 
is antisymmetric. However, the processes and -W^ are equally 

likely, so YQ(t")VQ(t') = 0. It now follows that ^/^(t' )Vj(t") = Vq(t ) V^(t ) 

exp(- I t' - t" 1), where t = min(t' ,t") . But it also can now be seen that, 
when t ' = t" = T , 


d V^CPVq(TT 
dt 


Vi(t)Vq(t) 
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Since Vj(0) = 0, V^(t ) Vq(t ) = 0 and consequently VQ(t‘ )Vj^(t") = 0 - as 
was to be shown. 

The asymptotic moments for which the lowest order terms were given are 

Y = + O(e^) 


r = ^ + O(e^) 

Note, however, that to this order Y is not Gaussian distributed. This is 
the case because the right side of equation (IV-23) is a nonlinear function 
of Gaussian random variables. 


Asymptotic Behavior of Langevin Model for Nonhomogeneous Turbulence 

In studies of turbulent dispersion the eddy diffusion model has repeat- 
edly been shown to be a good model where it is valid; that is, when in some 
sense t » T[^. The two comparisons between eddy diffusion and data for 
wall sources in chapter III support this observation. However, eddy diffu- 
sion generally is not valid at times that are short compared with Ti_; fig- 
ure IV-2 illustrates this. On the other hand, the discussion leading to the 
Langevin model was based only on an assumption of high turbulence Reynolds 
number; no restrictions were made on time. Therefore the Langevin model is 
uniformly valid in time. 

Because the eddy diffusion model agrees well with those experimental 
data to which it is applicable, our confidence in the Langevin model will be 
bolstered if it can be shown that this model reduces asymptotically to eddy 
diffusion when t >> T(_. If this can be shown, then eddy diffusion will 
be seen to be a special case of a high-Reynolds-number model. The Langevin 
model also will then seem to be the natural extension of eddy diffusion to 
the many situations where the latter is unsuitable. 

For the purpose of examining the asymptotic behavior of the Langevin 
model of dispersion in nonhomogeneous turbulence, we again consider the case 
of a slowly varying medium. Thus we consider the asymptotic behavior of 
equations (IV-22) and (IV-23) as t/T|_ . This asymptote can be 

obtained conveniently by letting T|_ 0. In the appendix to this chap- 

ter it is shown that the asymptotic behavior of equations (IV-22) and (IV-23) 
is obtained formally by equating the right sides of these equations to zero. 
By doing so and substituting dY = V dt, they become 

dYg = dW^ (IV-27a) 


^L 2 
= 2 Tl "^^0 ■ 




y/^i 


Yq dW, 


(IV-27b) 


61 



2 

Here it is necessary that YqVq dt = dYg/2 be substituted into equation 

(IV-23) before the asymptotic limit is taken. 

The solution to equation (IV-27b) subject to Yq(0) =0 is 


Vo - V^L »v“t 


(IV-28) 


2 2 

Hence Yg = 2T^^a^t, in agreement with the equation above equation (IV-10). 
When we substitute equation (IV-28), the equation for Y^ becomes 


dYl = 


pi p pi 


dW, 


where Ito's theorem 


“ 2 'l 


is used to write 


Yi = 



dt)/2. Hence 


(IV-29) 


— 2 ' 

Averaging gives Y^^ = OyT|_t, in agreement with equation (IV-26). Also 

T 

YgY^ = 0 (since =0), as was shown to be true generally for 

equations (IV-22) and (IV-23). 

Since Y = Yq + eY^ + ... , if equations (IV-28) and (IV-29) are 
added and the result differentiated, we find 

dY = dt + a^( + O(e^) 


But 



dt + 



(IV-30) 





eJ (0) 

- Yg + 

y2T[(0T ° 


0(€^) 


which is the bracketed expression in equation (IV-30). Hence to order 
equation (IV-30) can be written 


dY = a^TjY(t)]dt + yajTjY(t)]dW ^ (IV-31) 
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I 


This is just equation (III-l) for the eddy diffusivity 

K = o^lLCy) (IV-32) 

Thus we have achieved our objective of showing that the Langevin model 
(IV-19) tends asymptotically to the eddy diffusion model (III-l). It follows 
that the Langevin model is a plausible extension of eddy diffusion to disper- 
sion with finite time scale. Of course, our demonstration was limited to the 
case of slowly varying T[_ and was carried out only to order e; refer- 
ence 19 carries out the demonstration to order and also shows that the 
asymptote (IV-31) is generally correct. 


Case with Variable oy 

In equation (IV-19) ay was taken to be constant, but it was promised 
that later oy would be allowed to be a function of y. We now have a 
criterion for determining the form equation (IV-19) takes when oy is a 
function of y: the model (IV-19) must tend asymptotically to the diffusion 

model (IV-31) when T|_ + 0. 

It can be seen by following the procedure just used to derive equation 
(IV-31) that the form (IV-19) is not correct if ay is a function of y. 
However, since ay is constant there, that equation may be divided by 
2 

a^ and written 



2 

In this form a^T|^ plays the role that T^ played in equation (IV-19), and 

it can be shown that equation (IV-31) follows (ref. 19). Hence, for variable 
ay an appropriate Langevin model is^ 


^Since writing this report I have been told that Legg and Raupach (to appear in Boundary 
Layer Meteorology) use 


dV = - rr- dt + o 
\ 




to obtain the correct asymptotic behavior when is variable. Both this and equation 
(IV-33) are incorrect when the turbulence is homogeneous but nonstationary [tJy(t)J. 

The equation 







appears to be correct for nonhomogeneous and nonstationary turbulence. Clearly the issue is not 
finally resolved. 
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V 

"F 


(Y) 


V 


dt + 




dW, 


(IV-33) 


where K(y) = a^(y)\(y). 

An analysis of this model can be carried out for a slowly varying 
medium, as was done for the case with varying time scale. Define ')^(t) by 

y^(t) = (IV-34) 

ol CY(t)] 


and expand 

^(t) = r^it) + € r^{t) . . . 

% = °v^0) 2eay{0)a^(0)YQ . . . 


with Y and expanded as before in equations (IV-20) and (IV-21). Again 
Yq is the solution to the constant-coefficient Langevin model. Its variance 
is given by equation (IV-24) with T|^ and evaluated at y = 0. 

To order e we now have 


d -r 
~~St 


1 





(IV-35) 


for the mean value of (see equation following eq. (IV-24)). Because of 

the definition (IV-34), dY = aJ(Y) o^dt. Hence 


and 


dY 




2— *!v^ 

dt °v 1 0^ dt 


(IV-36) 


it being understood that a is evaiuateo ai y = u. me raccor a -jr . 
equation (IV-36) is given by equation (IV-25). However, an additional term 
associated with the inhomogeneity of appears in equation (IV-36). Upon 
substituting equations (IV-24) and (IV-25), equation (IV-36) can be integra- 
ted, subject to Yi(0) = 0. Thus it is found that 


is evaluated at y = 0. The factor 


1 n 


''l ’ - 1 


-t/T, 


9 ' 


-t/T, 


- 1 + ^ e 


-t/T, 


(IV-37) 
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By design, when t » Tj^ this reproduces equation (IV-26) with the eddy dif 
fusivity K = As t/T^^ 0, -*■ t^/2. Thus at short times the 

2 

mean drift is due to variation of and not to variation of Tj^; the vari 

3 3 

ation of T|^ enters at order t 11 ^. 

Perhaps this discussion of the Langevin model for nonhomogeneous dis- 
persion should conclude with the observation that this model has not been 
widely exploited. It seems a promising model for many problems involving 
dispersion with a finite time scale. 
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APPENDIX: ASYMPTOTE OF LAN6EVIN EQUATION 


Derivations of the asymptotic behavior of the Langevin equation are 
available in references 20 and 6. Essentially, a derivation is as follows: 
Since Tl is constant in equation (IV-1), it can be integrated as 


V(t) = V(0) 




(f-t)/TL 
e a 



(IV-Al) 


Because we are taking the limit t/T[_ ♦ ”, the first term can be dropped. 
Then 


Y 



V(t') dt' 



(t"-t')/TL 



(IV-A2) 


As T|_ + 0 


-|f-t"|/TL 


6(t' 


t") 


so equation (IV-A2) tends to 



or 



(IV-A3) 


Equation (IV-A3) is called the Smoluchowsky equation^. Equation (IV-22) is 

simpler and more rigorous derivation of Y(t) as t/T^ °° makes use of 
theorem 1* p. 8 of ref. 4. That theorem says that any continuous stochastic process X(t) that 

is independent of future events and has X(t) - X(t*) = 0 and [X(t) - X(t')]^ = t - t* (t >t') 

for arbitrary t,t' is a Wiener process. But if X(t) = Y(t)/a^ with Y(t) the solution 

to the Langevin equation, the properties of Y(t) derived in the text (eqs. (IV-9) and (IV-12) 
show that, as t/T^ X(t) has the required properties. 
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the same equation as equation (IV-1), so equation (IV-27a) follows irnmedi- 
ately, while equation (IV-23) can be integrated as 


V^(t) = 


■/' 


(t'-t)/T, 






(IV-A4) 


and equation (IV-27b) can be obtained by the argument leading to equation 
{IV-A3). Note that does not have a proper asymptote. (Its asymptote is 

not |/2Tj^ dW^/dt, as one might suggest from equation (IV-A3), for it was 

shown in the text after equation (IV-1) that Vq is statistically stationary 

with variance while (a^ ^2T|^ dW^/dt)^ = 2oyT|^/dt.) This is why it 

was necessary to write YqVq dt as dYQ/2 before taking the asymptotic 
limit of equation (IV-A4). 


67 



7 


REFERENCES 

1. Csanady, G. T.: Turbulent Diffusion in the Environment. D. Reidel 

Publishing Co., 1973. 

2. Tennekes, H.; and Lumley, J. L.: A First Course in Turbulence. MIT 

Press, 1972. 

3. Arnold, L.: Stochastic Differential Equations: Theory and Applica- 

tions. John Wiley & Sons, 1974. 

4. Gikhman, I. I.; and Skorohod, A. V.: Stochastic Differential Equations, 

Springer-Verlag, 1972. 

5. Monin, A. S.; and Yaglom, A. M.: Statistical Fluid Mechanics. Vol. I 

and II. MIT Press, 1975. 

6. Chandrasekhar, S.: Stochastic Problems in Physics and Astronomy. Rev. 

Mod. Phys., vol. 15, no. 1, Jan. 1943, pp. 1-89. 

7. Wax, N., ed.: Selected Papers on Noise and Stochastic Processes. Dover 

Publications, 1954. 

8. Mood, A. M.; Graybill, F. A.; and Boes, D. C.: Introduction to the 

Theory of Statistics. 3rd ed. McGraw-Hill, 1974. 

9. Rao, N. J.; Borwankar, J. D.; and Ramkrishna, D.: Numerical Solution 

of Ito Integral Equations, SIAM Journal of Control, vol. 12, no. 1, 

Feb. 1974, pp. 124-139. 

10. Mil'shtein, G. N.: Approximate Integration of Stochastic Differential 

Equations. Theory Probab. Its Appl . (Engl. Trans!.), vol. 19, 1974, 
pp. 557-562. 

11. Chatwin, P. C.: The Dispersion of a Puff of Passive Contaminant in the 

Constant Stress Region. Q. J. R. Meteorol . Soc., vol. 94, 1968, pp. 
350-360. 

12. Durbin, P. A.; and Hunt, J. C. R.: Dispersion from Elevated Sources in 

Turbulent Boundary Layers, J. Mech., vol. 19, no. 4, 1980, pp. 679-695. 

13. Shlien, D. J.; and Corrsin, S.: Dispersion Measurements in a Turbulent 

Boundary Layer. Int. J. Heat Mass Transfer, vol. 19, no. 4, Apr. 1976, 
pp. 285-295. 

14. Sullivan, P. J.: Longitudinal Dispersion Within a Two-Dimensional Tur- 

bulent Shear Flow. J. Fluid Mech., vol. 49, pt. 3, Oct. 1971, pp. 
551-576. 

15. Lighthill, M. J.: Introduction to Fourier Analysis and Generalised 

Functions. Cambridge University Press, 1975. 

16. Taylor, G. I.: The Present Position in the Theory of Turbulent Diffu- 

sion. Advances in Geophysics, vol. 6, H. E. Landsberg and J. VanMieghen, 
eds.. Academic Press, 1959, pp. 101-111. 


68 



17. Taylor, G. I.: Diffusion by Continuous Movements. London Mathematical 

Society Proceedings, vol. 20, ser. 2, 1921, pp. 196-212. 

18. Reid, 0. D.; Markov Chain Simulations of Vertical Dispersion in the 
Neutral Surface Layer for Surface and Elevated Releases. Boundary Layer 
Meteorol., vol. 16, no. 1, 1979, pp. 3-22. 

19. Durbin, P. A.: Asymptotic Behavior of the Langevin Equation in a Non- 

homogeneous Medium, unpublished (1982). 

20. Schuss, Z.: Singular Perturbation Methods in Stochastic Differential 

Equations of Mathematical Physics. SIAM Rev., vol. 22, no. 2, Apr. 

1980, pp. 119-155. 


69 



1. Report No. 

NASA RP-1103 


2. Government Accession No. 


3. Recipient's Catalog No. 


4. Title and Subtitle 

STOCHASTIC DIFFERENTIAL EQUATIONS AND TURBULENT 
DISPERSION 


5. Report Date 
April 1983 


6. Performing Organization Code 
505-32-02 


7. Author{s) 

Paul A. Durbin 


8. Performing Organization Report No. 
E-1425 


10. Work Unit No. 


9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and ^ace Administration 
Washington, D. C. 20546 


13. Type of Report and Period Covered 
Reference Publication 


14. Sponsoring Agency Code 


15. Supplementary Notes 


16. Abstract 

This report introduces aspects of the theory of continuous stochastic processes that seem to contribute to 
an understanding of turbulent dispersion. It is expected that the reader has a knowledge of basic probability 
theory and some exposure to turbulence theory and the phenomena of turbulent transport. The report is 
addressed to researchers in the subject of turbulent transport, especially those with an interest in stochastic 
modeling. However, the material presented herein deals with the theory and philosophy of modeling, rather 
than treating specific practical applications. The book by Csanady (ref. 1) is a good place to begin an ex- 
ploration of applications, as well as to obtain a background to the contents of the present report. 


17. Key Words (Suggested by Author(s) ) 
Stochastic processes 
Turbulence 
Transport 

Contaminant concentration 


18. Distribution Statement 

Unclassified - unlimited 
STAR Category 34 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

73 

A04 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1983 









